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ABSTRACT 
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O- ■ At least 10-15% of nearby sun-like stars have known Jupiter-mass planets. In 



contrast, very few planets are found in mature open and globular clusters such 



as the Hyades and 47 Tuc. We explore here the possibility that this dichotomy 
is due to the post-formation disruption of planetary systems associated with the 
stellar encounters in long-lived clusters. One supporting piece of evidence for 
q . this scenario is the discovery of freely floating low-mass objects in star forming 

+3 • regions. We use two independent numerical approaches, a hybrid Monte Carlo 

and a direct iV-body method, to simulate the impact of the encounters. We 
show that the results of numerical simulations are in reasonable agreement with 
analytical determinations in the adiabatic and impulsive limits. They indicate 
that distant stellar encounters generally do not significantly modify the com- 
pact and nearly circular orbits. However, moderately close stellar encounters, 
which are likely to occur in dense clusters, can excite planets' orbital eccentric- 
ity and induce dynamical instability in systems which are closely packed with 
multiple planets. The disruption of planetary systems occurs primarily through 
occasional nearly parabolic, non adiabatic encounters, though eccentricity of the 
planets evolves through repeated hyperbolic adiabatic encounters which accumu- 
late small-amplitude changes. The detached planets are generally retained by 
the potential of their host clusters as free floaters in young stellar clusters such 
as a Orionis. We compute effective cross sections for the dissolution of plane- 
tary systems and show that, for all initial eccentricities, dissolution occurs on 
time scales which are longer than the dispersion of small stellar associations, but 
shorter than the age of typical open and globular clusters. Although it is much 
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more difficult to disrupt short-period planets, close encounters can excite modest 
eccentricity among them, such that subsequent tidal dissipation leads to orbital 
decay, tidal inflation, and even disruption of the close-in planets. 

Subject headings: globular cluster, extra-solar planets - solar system: evolution 

1. Introduction 

The most successful observational method for planet searches has been the method of 
radial velocity surveys. Among over 1000 target stars, Jupiter- mass planets, with period up 
to a few years, have been found around more than 10 % of them (Cumming et al. 2007). This 
fraction is likely to increase with follow-up high-precision measurements of trends in the radial 
velocity curves which reveal suggestive signals for many additional planetary companions. 

This prolific success in the discovery of planets around these field stars is in sharp 
contrast to the rarity of planets around stars in both open and globular clusters. Despite 
several attempts to search for planets in the Hyades (Guenther et al. 2005), only one planet 
is found around one star (Sato et al. 2007). It is tempting to attribute this dichotomy to an 
environmental influence on the proficiency of planet formation. For example, Bonnell et al. 
(2001) suggested that, if the cluster went through an initial high density phase, close stellar 
encounters during the planetary formation epoch may have tidally truncated the disk and 
eliminated the domain of giant planet formation. This argument, however, is weakened by 
the poorly known time scale and location of giant planet formation. It also does not take 
into account the possibility of subsequent stellar accretion of the tidally stripped gas which 
may provide a protracted environment for planet formation. The discovery of planets in 
several binary systems with separation of a few AU's (Queloz et al. 2000, Hatzes et al. 2003, 
Zucker et al. 2004, Correia et al. 2007) is another sign that perturbations which are induced 
by neighboring stars alone cannot be an efficient mechanism to quench planet formation. 

The current paradigm for star formation theory is based on the realization that field 
stars are generally formed in dense gaseous environments which are similar to, but slightly 
less massive and concentrated than the progenitor of the longer-lived open clusters. In many 
dense star-forming molecular clouds, the fraction of stars with protostellar disks depends 
more sensitively on the stellar age rather than the dynamical properties of their host clusters 
(Hillenbrand 2005). Even in dense regions which contain massive stars, such as the Orion 
nebula, the silhouettes of protostellar disks are commonly found (O'Dell et al. 1993). Traces 
of radioactive decay products in the most primitive chondritic meteorites provide strong 
evidence that the formation of the solar system was preceeded by a nearby supernova event 
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(Cameron & Truran 1977). Since the massive progenitors of supernovae are only found in 
relatively dense molecular clouds, it is likely that the solar system actually formed in a stellar 
aggregate with ~ 2000 members (Adams & Laughlin 2001). These circumstantial pieces of 
evidence support the conjecture that planets may not form less prolifically in open clusters 
than the progenitors of the present-day field stars. 

An alternative scenario for the rarity of planet detection in open clusters may be at- 
tributed to post-formation dynamics. There is supporting circumstantial evidence for dy- 
namical disruption of planetary systems during the infancy of their host stars in dense 
molecular clouds. A population of "freely floating" planets, with 5-15 Jupiter masses, has 
been found in the young cluster a Orionis (Zapatero-Osorio et al. 2000, Lucas et al. 2001). 
If these free floaters indeed formed near some host stars, a reliable census of their frequency 
in a wide range of star forming regions would provide a calibration on the critical conditions 
for planetary system disruption. However, this task is particularly challenging because their 
detection and distinction from bound planets by microlensing observations are generally 
difficult (c.f. Han 2006). 

In contrast to open clusters, stellar associations which provided the birth place for the 
known planet-bearing field stars probably dispersed under the combined action of internal 
dynamical relaxation, stellar evolution, and the tidal perturbation of their host galaxies over a 
time scale ~ 10 8 yr. During this epoch, the planetary-retention probability of individual stars 
is determined by a competition between cluster disruption and planetary-system break up. 
The planetary-system dichotomy would be explained if stellar interactions and encounters 
affect the stability and dynamical evolution of their companion planetary systems on a time 
scale longer than the associations' dispersal time scale but shorter than the age of typical 
open clusters. 

On the theoretical side, Smith & Bonnell (2001) inferred that only a small fraction of 
planets would become detached from their host stars during the characteristic life span of 
young clusters and association. Their results are obtained using the equations of restricted 
three body motion to approximate encounters between planetary systems with nearly circular 
orbits and single stars. They claimed that the few detached planets would escape from open 
and young clusters because the planets' recoil speed is generally large compared with the 
clusters' velocity dispersion. These conclusions are modified if the perturbers are binary 
stars. This problem has been systematically addressed with a Monte Carlo approach by 
Laughlin & Adams (1998) in which they focused their attention on the effect of four body 
scattering (a binary star, and a planet). They averaged the stellar density over the cluster 
to obtain effective cross sections for the disruption of a Jupiter-like planetary orbit by all 
the binary stars. Also Davies & Sigurdsson (2001) included in their Monte Carlo study 
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encounters with binaries. 

In these early studies, a restricted range of planetary orbits was under consideration. For 
example, Laughlin and Adams focused their attention on planets which were initially on a 5 
AU circular orbit around the host star. However, the known extra solar planets have a semi- 
major axis ranging from 0.05 AU to several AU. Those with periods longer than a week have a 
nearly uniform eccentricity distribution up to unity (cf. Cumming et al, in prep.). Although 
planets are probably formed in their nascent disks beyond the snow line with nearly circular 
orbits (Ida & Lin 2004), their eccentricity may be excited shortly after their formation 
through their interaction with their nascent disks (Goldreich & Sari 2003, Creswell et al. 
2007) or with each other (Rasio & Ford 1996, Lin & Ida 1997). Dynamical instability and 
relaxation is particularly effective in exciting the planets' eccentricity especially in closely- 
packed multiple-planet systems (Papaloizou & Terquem 2002, Zhou et al 2007, Juric & 
Tremaine 2007, Ford & Rasio 2007). 

For most non-disruptive encounters, the stellar perturbation nonetheless modifies the 
eccentricities of the planets. Individual planets with modestly eccentric orbits may be more 
vulnerable to disruption by subsequent stellar perturbations. Repeated stellar interactions 
may be particularly damaging to the survival of planetary systems. Observational data also 
indicate that a significant fraction of known planets, if not most of them, have additional 
planetary siblings around their host stars (Fischer et al. 2001). Although the present-day 
separation in these known systems is generally wide and they are intrinsically stable, stellar 
perturbation in the past can provide an additional avenue for inducing their dynamical 
relaxation. The cumulative excitation of eccentricity, which results from both close and 
distant single and binary star encounters, may also lead to dynamical instabilities in closely 
packed multiple planetary systems, resulting in dissolution of a planetary system or in the 
merger of some planets with their host star. The destabilization of these multiple-planet 
systems generally occurs through secular and nonlinear interactions over many orbital periods 
(Zhou et al. 2007). 

Taking these physical processes into account, there is a need to investigate the effect of 
stellar encounters on planetary systems with more general initial condition than previously 
considered. Preliminary studies (see e.g. the different results obtained from direct iV-body 
simulations of Hurley & Shara 2002) indicate that the planetary retention efficiency may be 
strongly modified if the planets attain a modest eccentricity shortly after their birth. Ideally, 
the dynamical interaction between multiple planets and cluster stars needs to be investigated. 
In reality, the vast range of possible orbital configuration and mass distribution make any 
attempts to represent potential planetary systems futile. The orbital periods of compact 
planetary systems are much shorter than the crossing time of the cluster. On the time 
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scale of a large number of orbital periods, multiple planets interact with each other through 
cumulative low-amplitude secular perturbations. The accuracy requirement for calculating 
the long-term behaviour of multiple-planet systems is much higher than that needed for 
the reliable simulation of star clusters. For long-term solar system dynamics, symplectic 
methods, using a generalized leap-frog, like the widely used Wisdom-Holman symplectic 
mapping method (Wisdom & Holman 1991, see also review by Duncan & Quinn 1993), are 
the best suited integration methods. They do not show secular errors in energy and angular 
momentum. In principle, the symplectic mapping methods can be used to treat the planets, 
while the dynamics and perturbation induced by their host stars is computed with an N- 
body scheme. However, in their standard implementation, the symplectic mapping methods 
require a constant time step, which is not compatible with the central motivation of the 
Ahmad-Cohen neighbour scheme. Another more practical approach to strongly reducing the 
planets' secular errors in a iV-body scheme is to enforce a time-symmetric scheme by making 
the time steps reversible through an iteration (Hut, Makino, & McMillan 1995, Funato et 
al. 1996, Kokubo, Yoshinaga & Makino 1998). Such schemes have not yet been used for 
long-term secular evolution of planetary systems. 

Because of these complications, it would be immensely challenging to consider the ef- 
fects of secular interactions between closely-spaced multiple planets concurrently with the 
dynamical evolution of the cluster. As an alternative, we may aim to disentangle the differ- 
ent physical effects, and we begin here by focusing solely on the influence of gravitational 
encounters on an otherwise unperturbed planetary orbit. Having derived cross sections and 
time scales for these processes as functions of relevant parameters of the planetary system 
and its environment, we will be in a much better position to assess the conditions under 
which internal and external perturbations of planetary systems will couple, and those in 
which they will act on very different time scales. 

This analysis can also be applied to address the issue of a lack of planets in globular 
clusters. In a recent search for transit events among stars in the core region of an old, metal- 
poor globular cluster, 47 Tuc, no short-period planets have been found (Weldrake 2007 and 
references therein), even though 17 such objects were expected to be discovered (Gilliland et 
al. 2000). One possible cause for the absence of planets in globular clusters such as 47 Tuc 
may be the suppression of their formation due to lack of heavy elements in their progenitor 
clouds (Ida & Lin 2005a). However, we note that the mean metallicity deficiency (compared 
with the solar value) of 47 Tuc is smaller than the spread in the dust mass among disks 
around T Tauri stars (Beckwith 1998). The detection of another planet around a pulsar 
in globular cluster M4 confirms that although their formation may be suppressed it is not 
prohibited around metal-deficient stars. However, its stellar allegiance may have been altered 
through an exchange/ capture event during a close encounter between its original parent and 
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its present host star (Sigurdsson et al. 2003). 

Short-period (a few days) planets reside deeply in the potential of their host stars and 
can only be perturbed by a very close encounter between stars. Using an order of magnitude 
estimate for the encounter time scale between a single star and a planetary system, Bonnell 
et al. (2001) argued that, in globular clusters, planets with semi major axis greater than 0.3 
AU may be detached from their host stars. This estimate is based on an extrapolation from 
the outcomes of encounters between systems of three bodies of comparable mass, which may 
not be appropriate for the li miting case in which one of these bodies (the planet) is much 



less massive than the others (IFregeau et al.ll2006l ). 



In dense globular clusters, occasional close stellar encounters can induce eccentricity 
excitation even for planets with relatively short periods (Davies & Sigurdsson 2001). In 
multiple-planet systems, eccentricity excitation of long-period planets can also affect close- 
in planets through the secular interaction between them (cf Mardling & Lin 2002, Nagasawa 
& Lin 2005). In addition, tidal dissipation inside both short-period planets and their host 
stars can damp their eccentricity as indicated by the negligible eccentricity observed among 
all planets with period < 6 — 7 days (cf Dobbs-Dixon et al. 2004). The combined influence of 
these effects heats the planetary interior and may cause the planet to inflate (Bodenheimer 
et al. 2001). Very close to their host stars, tidal inflation may be sufficiently strong to 
disrupt planets (Gu et al. 2003). As the stellar relaxation process leads to an increase in the 
density of background stars, their perturbation on their planetary companions intensifies. 
This process may be particularly effective in eliminating planetary companions of relatively 
low-mass stars in a globular cluster environment such as 47 Tuc. 

In order to 1) verify this induced disruption scenario for the differences between plane- 
tary systems around stars inside and outside stellar clusters, 2) quantitatively determine the 
critical condition for planetary retention, and 3) assess the impact of stellar perturbations 
on close-in planets, we consider, in this paper, the dynamical evolution of planet-bearing 
stars in a range of cluster environments. In contrast to previous investigations, we focus on 
the diverse dynamical structure of planetary systems rather than the distribution of stellar 
kinematics, mass, and multiplicity. For simplicity of analysis, we only consider interactions 
between single stars and planetary systems. In realistic clusters where their fraction is large, 
binary stars may indeed provide more effective perturbers than single stars for the break-up 
of planetary systems, especially during the early epoch of cluster evolution when the stellar 
density is relatively high. Nevertheless, interactions between single stars are also important 
for the long-term survival of planetary systems in star clusters. 



In order to resolve some of these outstanding issues in both open and globular cluster 
environments, we carry out a series of numerical simulations which are designed to investi- 
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gate the influence of stellar encounters on the dynamics of extrasolar planets with a wide 
variety of orbital properties. In contrast to previous studies we use a special iV-body code 
(NBODY6++: Spurzem 1999), which is based on the NBODY6 code by Aarseth (1999a,b, 
2003) and can be used on massively parallel computers. In addition we use the hybrid Monte 
Carlo model (HMC) of Giersz & Spurzem (2003) as an approximate model of star clusters 
with large particle numbers and many planetary systems. 

These simulations, though more time consuming than plain statistical approaches in 
which only encounters are modeled, have the advantage that they can take into account 
both the spatial and time variation of the stellar background in young, open, and globular 
clusters. In a previous series of experiments, Davies & Sigurdsson (2001) have studied a 
Monte Carlo model of planetary systems interacting with stars and binaries in a dense 
star cluster such as 47 Tuc. However they did not include the dynamical evolution of the 
cluster and just studied isolated three- and four-body encounters. (Also their coverage of 
parameters was smaller than ours.) Both the NBODY6++ and HMC code are well suited 
to follow a wide range of encounters and to overcome some of these shortcomings. 

The NBODY6++ scheme is also ideally suited to the inclusion of the dominant dynam- 
ical influence of binary stars and occasional hierarchical triple systems, but in this paper we 
only focus on interactions of single stars with planetary systems. Another advantage is the 
possibility of simulating the consequence of several successive stellar encounters. After the 
planetary eccentricities are slightly excited by the first encounters, their rate of increase may 
be accelerated during the subsequent encounters. 

In this paper, we first focus on a comparison between empirical cross sections for orbital 
changes and analytical estimates (using several approximations). This study is intended to 
clarify the extent to which the existing analytical cross sections can be used, and we also are 
interested to see from the numerical models where they cannot be used any more, e.g. for the 
case of liberation of planets, either by successions of weak encounters or single strong ones. 
Our main interest is to understand the physical mechanisms, and an improvement for more 
realistic environments (binaries, multi-planetary systems and their internal interactions) is 
the subject of future work. Section 2 summarizes some of the existing knowledge on analytical 
cross sections for changes of eccentricity and semi-major axis of a planetary system due to 
encounters. In Section 3 we describe the setup of the numerical simulations and Section 4 
contains a description and analysis of the results. We pay special attention to the changes 
of planetary orbits in a statistical way by searching for encounter events which affect their 
eccentricities and semi- major axes (see for comparison also Theuns 1996). Once identified, 
data on these events are accumulated, and the results are then binned and presented as 
empirical differential cross sections for the outcomes of encounters between single stars and 
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planetary systems; here we use methods developed for the Monte-Carlo models and those of 
Giersz & Spurzem (2003). Finally, we summarize the findings of these numerical experiments 
and discuss their astrophysical implications in Section 5. 

2. Analytical Results on Encounters with Planetary Systems 

2.1. General Analytical Approach 

We aim to study the statistics of encounters between a passing star (the third body) 
and a simple planetary system consisting of a single planet and its star. We focus on the 
change in eccentricity 5e and the relative change in binding energy A = 5e/e. These depend 
on all the initial parameters of the encounter, but the essential parameters are the distance 
of closest approach, r p , and the speed of the third body at infinity, V (see Figs. [CJandE]). 
If r p 3> a, where a is the semi- major axis of the planet, then the encounter is tidal. If V 
is much larger than the orbital speed of the planet then the encounter is impulsive, unless 
r p is so large that the time scale of the encounter becomes comparable with the period of 
the planet. At still larger r p the encounter is adiabatic, which means here that the angular 
speed of the planet is much larger than that of the passing star. If V is small compared 
to the orbital speed of the planet, then a tidal encounter is always adiabatic. When r p is 
very large the path of the passing star is hyperbolic with high eccentricity, but at smaller 
values (though still in the tidal regime r p ^> a) it is nearly parabolic. These distinctions are 
essential in the analytical interpretation of numerical data. 

Appendix A collects a number of formulae for A and 5e which are approximately valid 
deep within some of these regimes. In adiabatic regimes an important role is played by 
parameters e' and K, which serve to quantify these ideas. The first parameter is just the 
eccentricity of the relative orbit of the passing star. The parameter K is defined by 

K= iw*ny {1) 

V M n3 W { ' 

where r p denotes the distance of pericentre for the passing star, the mass of the binary 
is M 12 = m% + m 2 , the mass of the third star is m 3 , and M i23 = M 12 + m 3 . In our 
simulations, mi = m 3 and m 2 <C mi, and so we have K = (r p /a) 3//2 . K measures the ratio 
of the time scales involved, and a necessary condition for an adiabatic encounter in the sense 
explained above is K ^> 1. However, this criterion by itself would be sufficient only if the 
encounter with the third star were parabolic. If the encounter is hyperbolic the condition 
for an adiabatic encounter is given by Kj\fd + 1 ^> 1 (Heggie & Rasio 1996), in which the 
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hyperbolic eccentricity 

/ / -nT/2 \ 2 

(2) 



/ I ( pV 2 \ 2 



where p, V denote the impact parameter and relative velocity at infinity. With these pa- 
rameters we get for the minimum distance r p of the third star relative to the centre of mass 
of the planetary system 

The impact parameter p and distance of closest approach r p are related by 



Now we consider the orbital response of a planet as a consequence of stellar encounters. 
Due to perturbations by passing stars, the semi-major axis and eccentricity of a planet 
experiences successive changes 8a and 5e which correspond to changes in its binding energy 
and angular momentum per unit mass, according to the equations 

A = f = -ir < 5 > 

hi a 
5 J 1 5a eSe 

-J=2^-— 2 - (6) 
The orbital energy per unit mass E is inversely proportional to the semi-major axis, and for 
a circular orbit it is just minus half the squared velocity. Accompanying the changes in e 
and a are changes in i, u, and Q, but we pay no attention to these in this paper. 

The formulae in Appendix A give approximate expressions for A and Se for a single 
encounter, and depend on all the initial conditions. For statistical purposes it is more useful 
to consider instead cross sections, a. For example cr(A), for A > 0, would be defined as the 
cross section for encounters in which the relative change of binding energy exceeds A. This 
would be a function of V, a, the masses and (possibly) e, though usually we will consider 
cross sections averaged over the distribution of e. Often also we consider differential cross 
sections such as da/dA. 

Cross sections essentially involve the impact parameter p, while the size of the changes 
A and 5e is more directly related to the distance of closest approach, r p . There are two 
limiting cases: (i) strongly hyperbolic encounters without significant gravitational focusing, 
where V 2 3> GM^/p, so that e' 3> 1, r p 3> GM^/V 2 , and hence p oc r p ; and (ii) the 
case in which V 2 GM 12 s/p, e' w 1, and so r p <^ GM^/V 2 . In this second case we have 
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strong gravitational focusing and the second term in equation (j3J) above dominates, so that 
p oc ^Jfp. When V = the encounter is parabolic and e' — 1. 

A given planetary system may experience encounters of either type, but which type 
dominates depends essentially on V and a. We have performed simulations for two ranges of 
semi-major axes of planetary orbits, which we denote for reference as "hard" (0.03 to 5 AU) 
and "soft" (3-50 AU), cf. Table 1. As is shown in Figs. 1 and 2 our "soft" planets nearly 
exclusively experience encounters where the velocity V of the encountering star at infinity is 
larger than the orbital velocity of the planet t> rb- In the case of "hard" planets (in our sense) 
we find a large amount of encounters (but not all) with the opposite case V < f or b- Note 
th at our use of t h e term "hard" and "soft" for planetary systems differs from the definition 
of iFregeau et all (120061 ). They use the velocity v c as a threshold, such that for V ~ v c the 
binding energy of the planetary system is comparable to the kinetic energy of the incoming 
star. Such definition is equivalent to what is used for binaries in star clusters (Heggie 1975), 
where hard binaries have a clear tendency to become harder in close encounters. However, 
in our models, the planetary mass is very small compared to the stellar mass (777.2 *C mi), 
so we have v c ~ f or b^2/^i "C w or b- An inspection of our Figs. 1 and 2 shows that there 
are practically no encounters with V < v c . Therefore we choose our definition of "hard" 
and "soft" planets, which we will use henceforth without quotation marks, just motivated to 
distinguish two different initial sets of semi-major axes of planetary systems. Indeed, in our 
models there are effectively only two regimes. For our soft planets, V typically exceeds the 
orbital speed of the planet V > t> or b- These encounters are impulsive if r p is not too large. 
They become adiabatic and hyperbolic if r p is large enough. On the other hand, many (but 
not all) encounters of our .hard planets fall into the category V < t> or b- Here tidal encounters 
at small enough r p are parabolic, but hyperbolic encounters occur at large enough r p . 



2.2. Cross Sections for the Change of Eccentricity 

2.2.1. Impulsive encounters 

In the case of soft planets, the most important tidal encounters may be treated as 
impulsive. First we provide an estimate for the form of the cross section, and then state an 
accurate result. 

Let r be the time scale for the encounter, i.e. r = r p /V p , where V p is the velocity of the 
passing star at pericentre. Our assumption here is that r < t orb , where t or b is the orbital time 
of the planetary system, i.e. we are far from the fully adiabatic limit. Then we approximate 
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the change of velocity of the orbiting planet due to the perturbation by the passing star as 

. . AGm 3 a 

5v^tSo,= t (7) 
v v r p 

where 5a is the difference between the acceleration of the planet (by the perturber) and that 
of the parent star (by the perturber), which we estimate by a tidal approximation. 

Now the relative changes in the binding energy, s, and angular momentum, J, of the 
planetary system may be estimated by 8v/v, and so by equation flBj) this may also be used 
to estimate D, which we define as the change in the square of the eccentricity, D = S(e 2 ). 
(This is more convenient than Se itself in some later analysis.) So we have 



Sv I Gm 



v V aV p 2 \r p 

For the last step above we have used m 3 m ni\. We also use V 2 = V 2 + 2G(mi + m^)/r p . 
Thus 

The total cross section is a = 7rp 2 ~ nr 2 oc D" 1 , and so we get for the differential cross 



da 
dD 



oc D- 2 (10) 



For a proper estimate we may adopt the approach which was taken by Heggie (1975) 
for the computation of the cross section for the relative energy change. As shown in Ap- 
pendix IB.lt the result is that 

da _ (Gma 3 ) 1 / 2 
dD~ 2Cl ^Di 

where we have specialized to the case in which the two stars (the passing star and the sun 
in the planetary system) have the same mass m, and C\ is a constant which is defined in 
terms of a certain average, and evaluates numerically to C\ = 0.883 approximately. 



2C 1 ^" bU , (11) 

t t r-)2 ' V ' 



2.2.2. Adiabatic encounters 

Now we consider adiabatic tidal encounters. Using a first-order expansion, Heggie (1975) 
and Heggie & Rasio (1996) obtained the net secular (i.e. long-term relative to the binary's 
orbital period) changes of eccentricity Se for parabolic and hyperbolic encounters between a 
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single field star and a close binary star. We adapt their formula to our case, in which the 
masses m,\ of the planet's host star and m 3 of the approaching star are large compared to 
the planetary mass m 2 ; so we have M 12 = rri\ + m 2 ~ mi and M 12 3 ~ mi + 777,3. equation 
(7) of Heggie & Rasio (1996) yields for this case 



Se = — - 



777q 



-(- 

4 \mxMi 



1/2 



3/2 



g(e',n,u,i), 



'123/ VW (1 + e') 3 / 2 

where g is the function in curly brackets in eq.( 1A3l) in the present paper. It reduces to g 
tt sin(2f2) sin 2 i for parabolic encounters, so to order of magnitude we use g 
For highly hyperbolic encounters we have similarly g 
and we use e' oc p 
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1 in this case. 



e', f 3 (e') = e'/(l + e'f' 2 



,/-0.5 



1 



~ v r p oc K 2 ^. It follows that 



5e 
Se 



par 



hyp 



'OC 
(OC 



K~ l oc r" 3 / 2 



^-4/3 



oc r, 



(13) 



From equations (JT3D we can deduce approximate total cross sections for both limits, using 
er oc p 2 oc r 2 for the hyperbolic case, and a oc p 2 oc r p for the parabolic case. It follows that 



a 



par 



'OC 
'OC 



(Se) 
(Se) 



-2/3 
-1 



cr hyp oc [Oe) * (14) 

and so we recover in the parabolic case the form of the result in Heggie & Rasio (1996) (see 
below). From this, differential cross sections are computed by 

da 



d(Se) 
da 



d(Se) 



par 



hyp 



'OC 



(OC 



(Se 
(Se) 



,"5/3 



(15) 



The details are given in Appendix B, where the following accurate results are obtained. 
(Note, however, that these are differential cross sections for the change in e 2 , i.e. D = S(e 2 ).) 
In the case of extremely hyperbolic encounters with hard planets, the calculation is carried 
out in Appendix IB.2| and yields a formula of the same form as for impulsive encounters: 



da 
dD 



IQC2VG 



(16) 



3 VD 2 

where C2 = 0.5932 approximately. For parabol ic encounters, h o wever , the corresponding 
result was essentially given by equation (19) of iHeggie fc Rasiol (119961 ). As shown in Ap- 
pendix IB. 31 this leads to the differential cross section 

1 2 



^ = l(i 57r ) 2 / 3 
dD 2V 1 



n5 

6 



Gma 



(D) 



-5/3 



(17) 



The form of this cross section is different, because of gravitational focusing. 
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2.3. Cross sections for the change in energy 

2.3.1. Impulsive encounters 
We consider now the relative change in energy, A, and begin with i mpulsiv e tidal en- 



counters. As already mentioned, formulae were given by Section 4.2 of iHeggid (119751 ) for 
the change in energy in an impulsive tidal encounter, and all that is needed is to specialize 
to the case m\ = = m and ni2 — > 0, and to harmonize the notation. He gives (his 
equation (4.24)) a cross section for the change y in binding energy x. In the language of the 
present paper his cross section is our differential cross section da/dA. With A = 5e/e we 
readily find that 

da 2n y/ Gma 3 

dA = Y VA? ' ( } 

though Heggie used units in which G = 1 and so we have restored it in its proper place here. 

This formula makes no distinction of the sign of A: it app e ars th at hardening and 



softening encounters are equally likely. But Section 4.3 of IHeggid (119751 ) also showed that 
an application of the principle of detailed balance could be used to estimate the difference 
in the rate of hardening and softening encounters. A similar principle, ho wever, applies 



to diff erential cross sections, and we employ it here. From equation (A3) in iHeggie fc Hut 



(119931 ). which is expressed in slightly different notation, it is easy to see that 



^(A, a, V)a^V* = g(A', a', V')a'^V'\ (19) 

where we have made explicit the dependence of the cross section on the initial semi-major 
axis of the planet and the initial speed of the passing star. On the left side we are here 
considering encounters which begin with values a, V and end with values a' = a/(l + A) 
and V. The right side considers information on the time-reversal of such encounters, and so 
(1 + A)(l + A') = 1. In the present case, with 1712 = 0, we have V = V. Substituting for a' 
and V we deduce that 

^(W) = (l + Ar^(__|_, a/(1 + A) , n (20) 



Substituting from equation ([18]) we deduce that 



^(A, a, V) = (l + A)- 3 ^("A, a, V). (21) 



Finally, expanding for small |A|, we find that 



(22) 
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This interesting result shows that hardening encounters are rarer than softening ones. 
The same qualitative conclusion results from considering non-tidal impulsive encounters, 
i.e. impulsive encounters in which the passing star approaches the planet or its sun to a 
distance less than a. For comple teness we present here two results for this case drawn from 
equation (4.12) in iHeggid (119751 ) . but specialized to the present masses and translated into 
the notation of the present paper. They are 



da 
dA 



V 2 A 2 



3A J 



3AJ 



;i + a) 



-5/2 



A < 



A > 0. 



(23) 



Again we find that softening encounters have a higher cross section than hardening encoun- 
ters. 

Though these cross sections apply to impulsive encount e rs, wh ich is a regime of relevance 
for soft planets, it was pointed out by Section 5.1 of lHeggiel (119751 ) that the results should be 
roughly applicable also for non-tidal encounters with hard systems, except for a correction 
due to gravitational focusing. The passing star is accelerated by the star of the planetary 
system until its relative velocity becomes compar able with that of the planet. In the case of 
stellar-mass binaries, which was the application in lHeggid (119751 ). the outcome for hardening 
encounters (A > 0) is complicated by the possible capture of the passing star. This does 
not occur in the case of a planetary system in the limit m 2 — > 0, and so equation (1231 
does apply roughly (with th e aforemention e d corr ection) to planetary systems. It also helps 
to explain the discovery by iFregeau et al.l (120061 ) that a planetary system of intermediate 
binding energy softens on average (provided that it is not so hard that capture of the passing 
star is energetically possible). Intermediate refers here to a planetary sys t em in the regime 
v c < V < v orb, where v c is the critical velocity according to IFregeau et al.l (120061 ) and v or b is 
the orbital velocity of the planet. 



2.3.2. Adiabatic encounters 

We can discuss changes of the binding energy in much the same way that we discussed 
change s in ecc entricity in Section l2.2.2l For the parabolic case, expressions for A are given by 



Heggid (119751 ) and, more explicitly, by Roy & Haddow (2003), while Heggie (2005) extends 
these to the hyperbolic case. Again we here give only enough to explain the form of the cross 
sections, and defer further detail to the Appendix. 

Let e denote the binding energy of a planetary system, with semi-major axis a = 
Gmim2/2\e\, and m 2 <C mi. Then the relative binding energy change A = Se/e for an 
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encounter becomes 

A = -V^fi(e')# 1/2 exp {-^Kf 2 {e'Yj F(e,u,Q,nt ) (24) 

The definition of F(e,u,Q,nto) can be obtained by comparison with equation (1A9I) . This 
factor is obtained from the one given by Roy & Haddow (2003), equation (19) and Heggie 
(2005), equation (11) by taking out the factor a 2 . Thus our function F is of order unity, 
and depends only on the eccentricity, orientation and phase of the planetary orbit. For the 
purpose of our paper only an average over all possible values is of interest, and to order of 
magnitude we use here Fwl. Now we quote the two functions 



, 3 yV 2 - 1 - arccos(l/e') 

/2(e) = 27! V^W • (25) 

We discuss the functions for the parabolic case e' — 1 and for the extremely hyperbolic case 
e' ^> 1; for the first one we have /i(e') = ^(e') = 1, and in this case we reproduce the result 
of Roy & Haddow (2003). In the hyperbolic limit we have asymptotically fi(e') ~ e /_5//4 and 
f 2 (e') ~ e'^ 1 / 2 , and so we get 

A par « -V^K l ^exp(~K) 

A hyp « -v^e'-^K^exp^-^e'- 1 /^. (26 ) 

As before we have e' oc p ~ r p for the hyperbolic case, and in both cases K oc rp^ 2 by 
definition, whence it follows that 



A par oc -0Frf exp(-^rf) 



A hyp oc - V7rr p 1/2 exp f _ ^ r pj ■ ( 27 ) 

To compute differential cross sections we need to invert the function A(r p ), which can be 
done easily only if the exponential function dominates. In that case we have 

Cpar oc p 2 oc r p oc (In A) 2 ^ 3 

a hyp oc p 2 oc r 2 oc (In A) 2 , (28) 
where we think of A as positive. Thus differential cross sections are obtained in the forms 

^ a x( lnA )" /3 

GtA par A 
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da 
dA 



oc — In A. 

hyp A 



(29) 



It is shown in Appendix [C] that, to leading order, the full results are 

2-nGma 1 



da 

dA 
da 

dA 



par 



hyp 



V 2 A 
Tra 3 V 2 1 



(In A) 



-1/3 



Gm A 



In A. 



(30) 



3. Computational Methods and Initial Model Parameters 

We use our hybrid Monte Carlo (HMC) scheme as well as direct iV-body simulations 
to study the changes in the orbital elements of planetary systems induced by encounters in 
stellar clusters. The main limitations of the models are that all stars are single (except for 
the planetary companions), all have equal mass, and there is no stellar evolution. The main 
loss of generality is the choice of a specific range of semi-major axis for the planetary orbits. 
Though we argue below that the range adopted is astrophysically justified, it is one aim 
of the analytical work to show how the results might be generalized to values outside this 
range. We collect data for encounters between stars and planetary systems in different ways 
in the two models, as described in the following sections. The initial parameters of iV-body 
and HMC models are summarized in Table [TJ 



3.1. Direct iV-Body models 

The N-body code we use is the publicly available NBODY6++ version. The equations 
of motion of each planetary system can be treated by regularization. Encounters and per- 
turbations by flybys affect the orbital elements of these regularized pairs, and a sufficiently 
strong interaction with a passing star can dynamically dissolve the planetary systems. The 
only special refinement to mention is that we utilize the newly applied classical method with 
which secular errors in the integration of close binaries in stellar systems can be strongly 
reduced (Mikkola & Aarseth 1998). 

For all models, we adopt an isotropic Plummer model for the stars' initial phase space 
distribution function. This model provides a reasonable approximation for open cluster 
potentials. All models are in dynamical equilibrium initially. Our model units are such 
that G = 1, M = 1, E = —0.25, for the gravitational constant, i nitial total mass and 



energy, respectively (standard iV-body units. iHeggie fc Mathieul (119861 )). For all models, the 
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individual mass of stars thus scales with 1/N. Planetary masses are set to the constant 
value of 10~ 10 in iV-body units for both the iV-body and HMC models. This choice ensures 
that all our planets are practically massless compared to the stellar masses. Physical units 
are obtained by (a) assigning an individual stellar mass (m*) in units of solar masses (see 
Table [I]) to the stars and (b) defining 1 pc as one iV-body unit (this is the virial radius 
R = GM/4\E\). The length scaling law also fixes what is one iV-body unit in AU (2.27 • 10 5 
for a system with solar- mass stars and 1 pc identified with one iV-body unit). 

For all runs done here all stars have equal mass, and there are no stellar binaries initially. 
We place one and only one planet around N p of the stars. For the direct iV-body models we 
use 20,000 objects, N p = 1000 of which are planets, which are randomly attached to 19,000 
equal mass single stars, to form a planetary system. This is a smaller particle number than 
in the HMC runs (see below). Other properties, such as the initial Plummer model, and the 
absence of tidal fields, are exactly as in the case of Monte Carlo models. 

In all of our iV-body runs except two (runs E, EH) all planets have the same initial 
eccentricity e , and we do different runs with differing values of e (Models 1-6: see Tabled]). 
The eccentricities of our models range from 0.01 (in model 1) to 0.99 (in model 6). The 
use of several different models with different eccentricities has two effects: (i) it improves 
the statistical data for the relatively small particle number (and relatively small number of 
scattering events) in the direct iV-body simulations; (ii) it allows us to consider the depen- 
dence of scattering statistics on the planet's orbital eccentricity, and to facilitate a more 
straightforward comparison with analytic models (Heggie & Rasio 1996, Roy & Haddow 
2003, Heggie 2005). Runs E and EH, however, cover all initial eccentricities using a thermal 
eccentricity distribution /(e) = 2e, which is also adopted in the HMC models. For small 
e's, this distribution approximates the Rayleigh distribution which is expected as a conse- 
quence of dynamical instability and relaxation (Zhou et al. 2007). In all calculations, the 
initial phase and orientation of the planets' orbits, including the direction of their angular 
momentum vector and their periastron, are randomized in the obvious senses. 

The semi major axes of the planets are chosen with a logarithmic distribution, i.e. a 
constant dN p J dloga, between 3-50 AU for run E and between 0.03 - 5 AU for run EH. 
The initial distribution of semi-major axes in the runs E and EH corresponds to the soft 
and hard planetary systems in the HMC runs, respectively. Note that we denote these 
planetary systems as soft and hard to distinguish between the two types of runs (with 
physically wider and closer planetary systems). We would like to remind th e reader again 



here t hat our definition of hard and soft differs from the dynamical definition of lFregeau et al. 



(120061 ). According to the dynamical definition a planetary system is hard if V < 0.75v c , 



where v c is the critical velocity of the three-body system. We choose the orbital velocity 
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f or b of the planetary system as parameter to distinguish encounters. If we use for the 
average V 2 in our models twice the cluster velocity dispersion, we get in iV-body units 
V 2 / v orb = ^(°au/2.27 • 10 5 ), where we have again taken one iV-body length unit to be lpc 
and (Iau is the planet's semi major axes in astronomical units. In Tabled] we give the value 
of Vyv or b for our initial planetary systems. 

In choosing these parameters we have been guided by theoretical and observational 
considerations. The range of 3-50 AU of semi major axis represents the location where gas 
giant planets are most likely to form (Ida & Lin 2005). Migration due to protoplanet-disk 
interaction may repopulate the regions interior to 3 AU (Lin, Bodenheimer & Richardson 
1996) and dynamical instabilities could eject planets beyond 50 AU (Lin & Ida 1997). But 
most of the gas giant planets may remain near the location of their formation. While gas 
giants are most likely to have formed with nearly circular orbits, dynamical instabilities 
could excite their eccentricities to the point of ejection. 

Up to now, one or a few Jupiter-mass planets, with periods less than a few years, have 
been found around less than 10% of the nearby solar- type stars (Marcy et al. 2005), though 
planets with longer periods are expected to be more common (Trilling et al. 1998, Armitage 
et al. 2002, Cumming et al. 2007). The total mass of planets is too small to significantly 
perturb the internal dynamics of any stellar cluster. Thus, for the simulations to be presented 
here, we include N p planets with very small masses of order 10 -10 , such that the mass ratio of 
planetary to average stellar mass is 1.9 • 10~ 6 , comparable to the mass of terrestrial planets. 
The total and individual mass of our planets is small enough, that they do not contribute 
any significant dynamical feedback to the cluster of N s stars. 

As for the stellar system itself, our aim has been to represent a typical rich young star 
cluster such as the Orion region. However, most stars are formed in binary and multiple 
systems. In stellar clusters, the presence of binary stars can significantly speed up the 
relaxation process due to their larger mass (Gao et al. 1991). They also strongly enhance 
the frequency of close three- and four-body encounters due to their larger cross section. This 
will also affect planetary systems (Laughlin & Adams 1998), and the influence of interactions 
with binary stars on planetary systems will be investigated in future work. 
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Table 1: Initial parameters of iV-body and HMC models a 



Model 


N* 


m*(M Q ) 


N p 


a(AU) 


e 




V/v orh b 


Type 


1 


1.9 


10 4 


1 


10 3 


3-50 


0.01 




0.49 - 2.05 


iV-body 


2 


1.9 


10 4 


1 


10 3 


3-50 


0.1 




0.49 - 2.05 


iV-body 


3 


1.9 


10 4 


1 


10 3 


3-50 


0.3 




0.49 - 2.05 


iV-body 


4 


1.9 


10 4 


1 


10 3 


3-50 


0.6 




0.49 - 2.05 


iV-body 


5 


1.9 


10 4 


1 


10 3 


3-50 


0.9 




0.49 - 2.05 


iV-body 


6 


1.9 


10 4 


1 


10 3 


3-50 


0.99 




0.49 - 2.05 


iV-body 


E 


1.9 


10 4 


1 


10 3 


3-50 


/(e) = 


2e 


0.49 - 2.05 


iV-body 


EH 


1.9 


10 4 


1 


10 3 


0.03-5 


f(e) = 


2e 


0.05 - 0.65 


iV-body 


Soft 


3.0 


10 5 


1 


3.0 • 10 4 


3-50 


f(e) = 


2e 


1.98 - 8.13 


HMC 


Hard 


3.0 


10 5 


1 


3.0 • 10 4 


0.03-5 


f(e) = 


2e 


0.20 - 2.57 


HMC 



a N* is the total initial number of stars and planetary systems, m*(M Q ) is the mass 
of one star in solar units, N p is the initial number of planetary systems, a (AU) is 
the initial range of the semi- major axes in AU, and e$ is the initial eccentricity. 
b we use V/v orh = ^ N(a A u/2.27 ■ 10 5 ) to obtain an approximate range of values, 
cf. Section 3.1. 



There is one important technical aspect of these A^-body runs that remains to be de- 
scribed. While direct A^-body models contain less intrinsic approximations than other sim- 
plified models, such as the HMC model, for our purposes there is a drawback, because it is 
very difficult to identify isolated two-body encounters in an A^-body model. In fact it has 
been discussed, whether a real A^-body system's relaxation process can be described by the 
standard model of uncorrelated small angle two-body encounters between individual stars 
(Theuns 1996). Generally, it is possible to identify an encounter by checking the minimum 
distance to the closest neighbour of any given particle, to get the r p and the velocity V p at 
closest distance. However, it is very difficult to determine the proper initial parameters of 
an encounter, because the scintillation and fluctuation of the A^-body potential perturbs any 
orbit even at moderate distances. Despite all these factors it is possible operationally to 
determine encounter event data, similar to that and to be compared with the HMC model. 
In fixed time intervals of one A^-body time unit (approximately one half-mass initial crossing 
time) we monitor orbital elements of all planetary systems (e, a). If any one of them has 
changed by more than 5 ■ 10~ 7 since the previous time, we assume an encounter took place. 
We measure 5e and A, and thus have a data bank of encounters for the A^-body system 
similarly to that for the HMC runs. From a theoretical point of view there may be uncer- 
tainties about this procedure, but we don't know any better alternatives, and judging from 
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the results it seems to be a reasonable operational procedure. The value of the hyperbolic 
eccentricity e' cannot be determined a priori; however from comparison with analytical mod- 
els (see figures below) one can see that the deduced values of e' lie in a completely reasonable 
range. 

3.2. Hybrid Monte Carlo method 

We use the hybrid Monte Carlo (HMC) method developed by Spurzem & Giersz (1996), 
Giersz & Spurzem (2000, 2003) to model the evolution of a star cluster with a large number 
of stars and planetary systems. The latter are considered as if they were a binary; binaries 
are treated with the Monte Carlo scheme to follow their relaxation with each other and 
with single stars in the cluster, while single stars are described by the anisotropic gaseous 
model based on the Fokker-Planck approximation (Louis & Spurzem 1991). Close encounters 
between planetary systems and a single star are followed as in the cited papers using a direct 
few-body integrator employing regularization methods (cf. e.g. Mikkola 1997). For HMC 
the planet mass in N-body units is very small, 6.33E — 12 and the mass ratio of planets to 
stars is 1.9E-6 also, as in the case of iV-body runs, and consistent with the notion that our 
planets have a mass comparable to terrestrial planets. 

We describe results of two sets of models, each with 300,000 single stars, 30,000 of which 
initially have a planet (see Table [1]). The semi- major axes of planetary systems introduce 
another independent scale into the otherwise scale-free iV-body system. We have chosen 
initially to keep the planetary scale a (semimajor axes) constant relative to the scale radius 
of the stellar system, independently of N. As a consequence the size of planetary systems 
measured in iV-body units is the same for HMC and iV-body models. Note that this means 
the squared orbital velocity of planets scales as 1/N, where N is the particle number. In 
terms of cluster velocity dispersion our planetary systems become weaker (softer) with larger 
N. We will discuss the scaling behaviour of encounters between planetary systems and single 
stars in more detail in sections 5.3 and 5.4. 

All parameter choices for the planetary systems are analogous to those of the A-body 
simulations, as explained before, and some more details are given in the table. For each set 
of parameters, two models with independently generated initial phase space distributions 
are adopted to boost the statistical significance of our results. The runs were continued for 
approximately five initial half-mass crossing times. This is a small interval for relaxation, 
but already enough to sample interactions between planetary systems and single stars ap- 
propriately. The time step of the HMC code is small enough to resolve these very accurately 
in time, because it is given by the gaseous model code for the single star component. Despite 
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of small time steps it is taken care that every binary is treated by the Monte Carlo procedure 
for relaxation with the correct rate and no over-relax ation occurs. A de t ailed description 
of this and the initial setup for encounters is given in iGiersz fc Spurzeml (120031 ). Here, we 
will only summarize the main points of the setup procedure. Whether a planetary system 
actually suffers from an encounter with a star is determined randomly for each time-step 
and planetary system. First, a maximum impact parameter is chosen to cover all physically 
interesting cases, p max = 2000 AU. From p max we determine an encounter rate A" using the 
simple prescription A" = nAv re \, where n, A = vrp^ ax , t> rel are the local planetary number 
density, the geometrical cross section corresponding to p max , and the actual relative velocity 
of the star and planetary system chosen for an encounter. In every time step we determine 
whether an encounter takes place by comparing the probability obtained from the encounter 
rate with a random number. After we have identified an encounter we pick an actual impact 
parameter p, from a random distribution in p 2 . Using this with energy and momentum con- 
servation we transform the encounter from the cluster frame into the interaction frame and 
start the actual direct integration of the encounter using packages from Aarseth's NBODY6 
code. After the termination of the direct integration of the encounter the new parameters 
of the planetary system are computed and the outcome of the encounter is identified: flyby, 
dissolution or exchange and transformed from the interaction frame into the cluster frame. 



4. Results of numerical experiments 

The three parts of this section summarize our numerical results and their interpretation. 
We begin with a discussion of the destruction of planetary systems in stellar encounters, and 
the escape of the resulting population of freely-floating planets. Then we give a largely 
empirical description of the results of the non-destructive encounters. This analysis paves 
the way for the final subsection, in which we interpret the results as directly as possible in 
terms of differential cross sections. 

We note here that many of our encounters lead to very small changes of the eccentricity 
or semi-major axis of the planetary system (see some of the following plots). The changes are 
often sufficiently small as to cast doubt on the significance of such a result . Numerical errors 
are introduced in the three-body integration for HMC and by the stochastic background noise 
of potential fluctuations which are present in the direct A^-body simulations. To distinguish 
broadly those encounters which might be regarded as unreliable from the remaining results, 
we define certain criteria, namely K < 80 and |A| > 5- 10~ 7 , which help to identify the most 
robust results. In some of the following interpretative plots we show results from the full set 
of encounters, and in other cases only a limited subset. The selection of the representative 
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Table 2: Summary of results of hybrid Monte Carlo (HMC) runs 



Model 


HMC soft 1 


2 


HMC hard 1 


2 




1995262 


2426598 


3406206 


3650990 


t TV— body 


15.28 


15.68 


20.67 


30.31 


^~cr — ^ / ^cr 


5.40 


5.54 


7.30 


10.71 


r rh = t/txhfl 


5.93E-03 


6.08E-03 


8.02E-03 


1.18E-02 


-Wpl-disa 


149 


162 


7 


9 


-^pl— diss— esc 


10 


10 


3 


3 


A^pi-fr 


139 


152 


4 


6 


^Vpl-ff/^rh 


23440 


25000 


498.1 


508.5 




25.74 


27.44 


0.548 


0.560 


^pl-ff.rh 


0.7813 


0.8333 


0.0166 


0.0170 


•£pl— ff,cr 


8.58E-04 


9.15E-04 


1.83E-05 


1.87E-05 


Scaled £ p i-ff,rh 


0.0673 


0.0718 


1.43E-03 


1.46E-03 



-^events is the number of interactions with planetary systems, tv-body is the time of 
termination of the simulation in iV-body units, t cr $ is the initial crossing time, t r h,o is 
the initial half-mass relaxation time, ./Vpi_di ss is the number of dissolved planetary systems, 
N pl —diss— esc is the number of planets escaped from the system after the dissolution of the 
planetary system, Npis is the number of "freely floating" planets. 

The quantities x p is tT h and x p i_fj )Cr denote the probability for one planet to become 
a free floater, and are just obtained from the previous two lines by dividing by the 
total number of planets. x p i_fl ]Cr can be directly compared with iV-body results below, 
compare Eq. [Ml x p i_g irn will scale with the two-body relaxation time and therefore as 
N 9 .kihNi)/ NihihN 9 ), with iVi = 300.000 (HMC run), N 2 = 19.000 (iV-body run), 
7 = 0.11, cf. iGiersz Heggie! (E994h 



data will be clearly stated in the respective paragraphs or figure captions. 



4.1. Dissolution of planetary systems 

First, we examine the overall statistics on planetary-system retention. In Table [2] we 
provide some basic data which are generated with the Monte Carlo (HMC) scheme. The 
simulations were stopped after a few million encounters had taken place between planets and 
passing stars. The actual evolutionary duration corresponds to some 5-10 initial half-mass 
crossing times. The table summarizes some interesting information on the dissolution of 
planetary systems and the creation of free floaters in the HMC models. For comparison, we 
also present in Table [3j the analogous data which are generated with the iV-body scheme. 



-23- 



Table 3: Summary of results of direct iV-body runs 



Model 


1 


2 


3 


4 


5 


6 


E 


EH 


-^events 


24338 


45071 


52060 


70151 


86951 


81655 


62022 


6263 


^N— body 


72.0 


130. 


153. 


170. 


181. 


166. 


134. 


200. 


7~er 


25.5 


46.0 


54.1 


60.1 


64.0 


58.7 


47.4 


70.7 


Trh 


0.32 


0.59 


0.69 


0.77 


0.81 


0.75 


0.60 


0.9 


-/Vpl-diss 


33 


66 


69 


89 


63 


63 


70 


2 


Np\— diss— esc 


3 


5 


4 


5 


8 


8 


3 





iVpi-ff 


30 


61 


65 


84 


55 


55 


67 


2 


N p is/r rh 


93.8 


103.4 


94.2 


109.1 


67.9 


73.3 


111.7 


2.22 


N p is/t ci 


1.18 


1.33 


1.20 


1.39 


0.85 


0.94 


1.41 


0.028 


•£pl— ff,cr 


1.18E-03 1.33E-03 


1.20E-03 


1.39E-03 


8.5E-04 


9.4E-04 


1.41E-03 


2.80E-05 


£pl-ff,rh 


0.0938 


0.103 


0.0942 


0.109 


0.0679 


0.0733 


0.112 


2.22E-03 


£pl-ff,cr 


1.35 


1.48 


1.35 


1.57 


0.98 


1.05 


1.61 


1.54 


Cpl-ff,rh 


1.33 


1.50 


1.35 


1.57 


0.96 


1.06 


1.59 


1.51 



All quantities have the same meaning as in Table [2j The last two lines give £, defined as the ratio of 
x obtained from the iV-body model divided by the average of the two corresponding HMC results, 

i-e- £pl-ff,cr = £pl-ff,cr,Nbody/£pl-ff,cr,HMC,av and £pl-ff,rh = X p l_ff , r h,Nbody/^pl-ff,rh,HMC,av,scaled • 



For all models, we provide, in the tables, the rates of free floater liberation, alterna- 
tively in terms of the cluster's crossing or relaxation time. In the models which represent a 
population of soft planets, we create about one free floater per crossing time and about 100 
per relaxation time (iV-body). In contrast to previous claims (Smith & Bonnell 2001), only 
very few planets escape during the whole simulation - one order of magnitude less than there 
are free floaters, even though the initial orbital velocity of typical planets is larger than the 
velocity dispersion of the cluster. 

In order to rule out the possibility that the retention of disrupted planets in the shallow 
potential of the host cluster may be due to an under-representation of close encounters by 
the HMC scheme, an analogous model is simulated with the iV-body scheme. In order to 
compare the rates of free floater and escaper creation between these two approaches, which 
were obtained for practical and technical reasons with different numbers of particles and 
planetary systems, we first have to apply an appropriate scaling factor. In our system of 
iV-body units (used here for both methods) the total mass is constant (unity) and individual 
stellar masses scale as 1/N. Rates per crossing time scale only by the number of planets, 
because in iV-body units the cross section, if gravitational focusing prevails, scales inversely 
proportional to iV and the stellar number density is proportional to N. However, in the 
comparisons between the rates per relaxation time we should additionally scale by the ratio 
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of half-mass relaxation times for the two models. The exact scaling expression is given in 
the caption of Table [2j 

The tables show that both the rates per crossing time per planet, and the (scaled) rates 
per relaxation time agree well, but the numbers differ, in the sense that the rate of free floater 
creation is larger by some 35 % in the iV-body system. However, the intrinsic variation of 
iV-body results for different eccentricities is as large as that, therefore we conclude that our 
HMC and TV-body results do agree reasonably with each other. The iV-body result is that 
approximately 100 dissolutions of planetary systems occur per relaxation time in this case, 
regardless of their initial eccentricity. Contrary to the standard expectation that highly 
eccentric planetary systems should be more prone to disruption (Hurley & Shara 2002), we 
find that the disruption rate is essentially insensitive to the planets' original eccentricity. In 
fact, there is a slight tendency for fewer disruptions, in the case where all planetary systems 
are initialized with e = 0.99. 

In order to understand this puzzling result, we have examined the encounter activity 
between planetary systems and stars in the iV-body models. We find that, for highly eccentric 
planetary systems, there is indeed more encounter activity (e.g. more KS terminations) 
by some 40% (as compared to models with modest eccentricities, e.g. 0.6). Nevertheless, 
these encounters do not necessarily lead to free-floaters. After a brief episode of temporary 
liberation (KS termination), some of the highly eccentric planets quickly become attached 
to intruding stars. The exchange of the host stars occurs as a consequence of three- or 
more-body interaction. The many-body effects may also contribute to the small residual 
differences in the number of disrupting planetary systems between the HMC and iV-body 
simulations because, in general, the many-body effects should increase the dissolution rate 
(as seen in the iV-body simulations). This difference vanishes for the simulations of stellar 
encounters with hard planetary systems because the stellar cluster and the planetary systems 
are dynamically well segregated and the many-body effect should play a lesser role. While 
this effect is potentially very interesting, it is not in the scope of this paper to study it in 
further detail. 



4.2. Semi major axis and eccentricity changes 

With the aid of analytic formulae, we now analyse the consequences of stellar encounters 
including the vast majority which did not lead to the dissolution of the planetary systems. 
The main objective in this section is to determine the cross sections for eccentricity and semi 
major axis changes. 
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4-2.1. Domains of encounter classes 

In the analysis of the numerical data, it is useful to make direct comparison with the 
analytic results in §2. Following the prescriptions in §2 (Heggie 2005), we first show the 
location of the limited set of encounters in the normalized relative velocity-impact parameter 
plane for the hard and soft planetary systems (Figs. [Hand [2]). In these figures, the ordinate 
is V I a/ GM^/a, while the abscissa is given by r p /a for each encounter. The parameter space 
is separated by three domains in accordance with the nature of the encounters. The vertical 
line (r p = a) separates very close interactions from tidal encounters. The line separating 
adiabatic from non-adiabatic encounters is defined by V/r p = v c /a, where v c = \J GMi 2 / a is 
the circular velocity of the planet. The line which separates near-parabolic from hyperbolic 
encounters is given by the condition e' = 2. 

Figure. [1] represents a model for the hard planetary systems. Only a small number 
of encounters (for small r p ) with hard planetary systems are non-adiabatic. There is a 
considerable number of near-parabolic encounters for r p > a. In contrast, the cloud of 
representative points shifts upward in the model for the soft planetary systems (Fig. |2J). 
This systematic difference is due to the planets' larger semi-major axes. Consequently, there 
are many non-adiabatic encounters with hyperbolic speeds and r p > a for the soft planetary 
systems. The number of near-parabolic encounters is negligible. Very small changes in 
|A|(< 5 • 10~ 7 ) may be affected by numerical errors. The omission of these potentially 
spurious points accounts for the absence of points in the lower right corner of the figures. 
The condition that K < 80 also truncates the distribution of data points on the right hand 
side. A corresponding set of figures cannot be constructed for the iV-body simulations due 
to the lack of complete dynamical information. In these simulations the magnitude of v for 
individual encounters is not well determined. 



4-2.2. Correlated changes of orbital elements 

We now analyze the correlations between the relative change in the planet's energy and 
the magnitude of 5e. Figures. [3] and H] illustrate the full set of encounters for the soft planetary 
systems in the HMC and iV-body models (run E), respectively. Here we can clearly identify 
a huge number of encounters with very small changes. The data points on the upper right 
hand corner of these figures represent planetary systems with modest changes in both the 
relative energy and eccentricity. These figures show that there is a correlation between the 
relative semi-major axis and eccentricity changes. 

This correlation is much weaker for the bulk of encounters with very small changes, 
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Fig. 1. — Location of encounters for hard planetary systems in the hybrid Monte Carlo 
model, plotted according to the scaled velocity at infinity and the minimum distance r p 
in units of the semi-major axis a of the planetary system. Solid lines indicate boundaries 
between encounters which are close or wide, adiabatic or non-adiabatic, hyperbolic or near- 
parabolic. All details of scaling and definition of the boundaries are given in the main text: 
note that they are analogous to those in Fig. 1 of Heggie 2005. 

which is one of the reasons we consider changes at this level may be due to spurious random 
noise. For the iV-body models, the uncorrelated data representing very small orbital changes 
may also be due to the unresolved potential fluctuations and multiple concurrent encounters 
(see discussion above on how we sample the encounter data in the iV-body model). In the 
HMC model, where the maximum impact parameter in the simulation was set to be 2000 
AU, even the weakest encounters are well defined. The uncorrelated measurement of energy 
and eccentricity changes represents the small limits of numerical accuracy in the few-body 
integration and in their initialisation. 

As in the A^-body data, however, for the data analysis we only use those encounters 
which obey |A|, |5e| < 5 • 10™ 7 to exclude numerically unreliable results (see discussion in 
Sec. I3.ip . Both figures show that the changes are very symmetric with respect to their sign, 
i.e. there are as many cases with eccentricity decrease as with eccentricity increase. Note 
that for the system with larger N (the HMC models have 300,000 stars, as compared to 
the A^-body models, which had only 19,000), the median values of A and 5e are roughly 
one order of magnitude smaller. This difference is consistent with the interpretation that 
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Fig. 2. — Same as Fig. [TJ but for soft planetary systems. 
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Fig. 3. — Relative energy change vs. eccentricity change for soft planetary systems in the 
hybrid Monte Carlo model. All encounters are included. 

encounters become weaker for larger systems. 
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Fig. 4. — As Fig. [3J but for the direct iV-body model (soft-planetary systems) with an 
initially thermal eccentricity distribution. 



4-2.3. Parabolic versus hyperbolic encounters 

The next three Figs. El El and [7] again illustrate the dynamical changes in the plane 
defined by A and 5e. The main purpose of this analysis is to show the dependence of the 
outcome on the initial relative speeds between the encountering stars. Indeed it was the 
study of these results which forced us to extend theoretical results in the literature to the 
case of hyperbolic encounters. 

For the HMC model (Fig. ED we show three sets of encounters, depicted by different 
colours of points: e' < 2 (red), 10 < e' < 30 (green), and ef > 50 (blue); the general 
limits defined above (K < 80, A > 5 x 10 -7 ) were also applied. These three groups of data 
points represent nearly parabolic, intermediate, and hyperbolic encounters. Note that most 
hyperbolic encounters lead to very small changes in a and e. Only the parabolic encounters 
satisfy the necessary condition for planetary disruption and major orbital element changes, 
i.e. |A| > 1 and be ~ 1 respectively. 

For the iV-body model (Figs. El and [7]), although the encounters cannot be individually 
distinguished by their e', the overall results agree with those obtained from the HMC model. 
(Here data for very small changes, i.e. with 5e < 5 x 10~ 6 , are omitted, see figure). Figure El 
shows iV-body results with an initially thermal eccentricity distribution, while Fig. [7| presents 
results from run 1 (with e® = 0.01 initially for all planetary systems). 
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Fig. 5. — Location of encounters for hard planetary systems in the hybrid Monte Carlo model, 
plotted by relative energy change and eccentricity change, for three cases with different e' 
(near-parabolic, intermediate, and extremely hyperbolic; see main text for more details). 
For comparison solid lines are plotted from analytic expressions of Heggie & Rasio (1996) 
and Heggie (2005). 

In order to compare these numerical data with the analytic results (see §2), we super 
impose in all three figures, the correlated magnitude of d~e and A derived from the analytic 
expressions in equations (lT2"j) and (|23J), for three different values of e' = 1, 20, 80, respectively. 
The results of the HMC simulations and those of the iV-body model E (with an initial thermal 
eccentricity distribution), i.e. Figs. [5] and [6], agree fairly well with each other and with theory. 
Both models indicate that planets break up primarily through parabolic encounters. Note 
that e' increases from right to left in all three figures. It can also be seen that the encounters 
which have been excluded (due to small changes or large K) are extremely hyperbolic. 

4-2.4- Planetary systems with nearly circular initial orbits 

This case requires separate discussion. The plot of the iV-body results for eo = 0.01 
(nearly circular orbits, Fig. [7j) shows some differences from the previous two cases. For 
planetary systems which suffer large orbital changes, the correlation between 5e and A is 
much more pronounced, with less scattering of the points at the upper right hand side of the 
figure. For a given d~e these encounters for nearly circular binaries tend to exhibit smaller 
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Fig. 6. — Location of encounters for planetary systems with an initially thermal eccentricity 
distribution in the direct iV-body model (soft-planetary systems), plotted by relative energy 
change and eccentricity change. For comparison the same solid lines as in Fig. [5] are plotted, 
based on analytic expressions of Heggie & Rasio (1996) and Heggie (2005). 

energy changes than in the previous cases. 

Now we turn to a theoretical interpretation of these differences. Several theoretical 
results change when one considers strictly circular initial orbits (cf. Appendix lA.il) . This can 
also be seen in Fig. [7J where we have plotted curves based on assuming the same dependence 
of 5e and A on e' as for non-circular orbits. This comparison is unsatisfactory, particularly 
for the small values of K, which correspond to the limit of strong encounters. 

First consider the encounters for which |<5e| ^> e = 0.01. In this case analytical theory 
for a circular binary is relevant. In eqs. (1A4I) and (1A5|) . which give the change in eccentricity 
for this case, the main dependence on the distance of closest approach is in the exponential 
factors. The same exponential dependence occurs in the formula for the relative energy 
change in the parabolic case (eq. (1A10I) ) and in the hyperbolic case (not shown). Therefore 
we may expect A oc be for large changes in eccentricity (and energy). The dashed line in 
FigJTJhas the same slope, and it is followed by the trend of the points for large Se and A. 

There is a further complication, because the figure we are discussing deals with soft 
planetary systems, and the closest encounters are distinctly non-adiabatic (Fig [5]), and indeed 
impulsive. To find out how 5e and A are correlated in this case we may make use of eqs. ([5]), 
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([6]), from which it can be deduced that 




SE S(J 2 ) 



E J 2 



Now for an impulsive encounter with a circular binary 5E = v.5v and 5 J 2 = 2a 2 v.<5v. But 
we also have J 2 = GMa in this case, where M is the total mass of the binary, and a quick 
calculation shows that the lowest order contribution to S(e 2 ) vanishes. (It is obvious that it 
must do, as this contribution would be proportional to v.Sv, and would give both positive 
and negative values, whereas we must have S(e 2 ) > for an initially circular binary.) It 
follows that S(e 2 ) must be proportional to |5v| 2 , in this case, and so Se oc \Sv\. But SE is 
also proportional to \Sv\, as mentioned in Sec l2.2.1[ and so again we deduce that A oc Se. 

In this case of impulsive encounters, it was mentioned in Sec J2.3.11 that A ~ D = S(e 2 ). 
This suggests a steepening of the trend in the dependence of A on Se at the extreme right 
of the diagram, and this is not observed. 

In the remainder of the diagram (\Se\ <C 0.01) the binary essentially behaves as one 
with a non-zero initial eccentricity, and the distribution of points may be expected to behave 
much as in the two previous cases. The analytical expressions for Se and A, i.e. eqs. (lA3p 
and (lA9p . are approximately proportional to e for small e, and so the curves in this diagram 
should be scaled by a factor of order 0.01. This brings them into better accord with the points 
corresponding to more distant encounters. The numerical results indicate that the amplitude 
of eccentricity change is considerably larger than that of the change in energy or semi major 
axis. Since the dissolution rate is essentially independent of the initial eccentricity, planets 
which formed with circular orbits can quickly become eccentric. 



In §2, we show that the magnitude of a planet's dynamical response to a stellar encounter 
depends on both the impact speed (through e') and the distance of closest approach (through 
K). We now interpret our numerical results in terms of the dependence of K given by the 
analytic expressions. In Figure [8] the dynamical changes in the planetary orbits are depicted 
as in the three figures before, but this time plotted in the plane of K and A. Only HMC 
simulations are represented, for A^-body models it was impossible to obtain reliable encounter 
parameters for each interaction and compute K (see Sec. 13.11) . This plot can be compared 
with the analytic results more directly because the latter are often expressed as functions of 
K, whereas the correlation with Se is more indirect. 



4-2.5. Dependence on the intruding star's distance of closest approach 



The numerical results obtained with the HMC scheme (Fig. [8]) agree very well with 
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Fig. 7. — Location of encounters for planetary systems with initial eccentricity eo = 0.01 
in the direct iV-body model for soft planetary systems, plotted by relative energy change 
and eccentricity change. Again the same solid lines as in Fig. [5] are plotted from analytic 
expressions of Heggie & Rasio (1996), Roy & Haddow (2003) and Heggie (2005). 



the analytic expression. In this case, the fractional energy change at large K is correlated 
with e'. This trend indicates that in the high K limit, a vast majority of encounters are 
very hyperbolic. Only nearly parabolic (e' < 2) and non adiabatic, close (small K) stellar 
encounters lead to planets' disruption, i.e. |A| > 1. In addition, it follows from Fig. [5] that 
most non adiabatic, close encounters lead to major eccentricity changes, i.e. \5e\ ~ 1. 



4.3. Rates of energy and eccentricity evolution 

4-3.1. Distributions of 5 e and A 

For another comparison of the HMC and NBODY6++ results we have just binned the 
changes of A and 5e and computed the numbers N(A), N(5e) normalized to the total number 
of events. Figs. 191 and [TUl show the results as a function of 5e and A, respectively. Note that 
this is an unnormalized measure of the total cross sections. In the following section we will 
deduce properly the normalized differential cross sections as well. 

Here, we first observe that the shapes of the encounter frequency distributions are 
similar, albeit they are shifted (the numbers for the HMC model are smaller by about one 
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Fig. 8. — Location of encounters for hard planetary systems in the hybrid Monte Carlo model, 
plotted by relative energy change as a function of K = (r p / a) 3//2 , for three cases with different 
e' (near-parabolic, intermediate, extremely hyperbolic; see main text for more details). For 
comparison solid lines based on analytic expressions of Heggie (2005) are plotted. 



order of magnitude, except for small \be\). There are two possible explanations for this 
shift, and both effects may contribute. With a larger particle number in the HMC model 
all changes are on average about an order of magnitude smaller than in the iV-body model 
(compare discussion of Fig. [2D and so the curves are shifted to the left. On the other 
hand, in the iV-body system individual encounters with extremely small changes (which are 
easily observed in the HMC model) are obscured by the background of stochastic potential 
fluctuations. 

The nearly logarithmic distribution N(5e) oc be^ 1 for be > 10~ 4 implies that the mean 
square changes of eccentricity are dominated by close encounters and large be. We also note 
here that, in the HMC models we observe that the initially thermal eccentricity distribution 
is preserved, as would be expected. 

Turning now to the distribution of A, we note that, for energy changes with A ~ 0.1, 
negative changes of A are more probable, i.e. that ther e is a preferred t rend t owards softening 



of planetary orbits. This confirms the discovery by iFregeau et al.l (120061 ) that there is a 
range of intermediate planetary orbits with v c < V < v ov \ } for which there is a net surplus of 
softening encounters as compared to hardening ones. We also confirm this with our analytical 
results, as discussed qualitatively at the end of Section 2.3.1, and the result is also visible in 
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Fig. 9. — Fractional number of encounters as a function of eccentricity change (soft planetary 
systems), compared between the hybrid Monte Carlo and iV-body models with an initially 
thermal eccentricity distribution. The data are given separately for positive and negative 
changes. 



Fig. 13 (differential cross sections) with better statistical quality. 



4-3.2. Normalization for differential cross sections 

Finally we compute properly normalized differential cross sections, and compare them 
with theory. We restrict the numerical data to those from the HMC runs. Differential cross 
sections can be obtained from our numerical results using the binned data N(A) described 
previously. Because the cross sections depend on a and V, however (Section [2]), some post- 
processing is necessary. 

We define P to be the probability that a given encounter results in a value of A which 
lies within a certain bin. Then 

1 da 



P=-—-(A,a,V)5A, (31) 

n Pmax 

where 5 A is the range of the bin. Let us suppose that the differential cross section has a 
power-law dependence on a and V, i.e. 



§(W) = ^(A,MK^ 



(32) 
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Fig. 10. — Fractional number of encounters as a function of relative energy change (soft 
planetary systems), compared between the hybrid Monte Carlo and iV-body models with 
an initially thermal eccentricity distribution. The data are given separately for positive and 
negative changes. 

where a,f3 are certain constants and da/dA(A, 1, 1) is the normalization factor. Then, 
equating N(A) to its expected value, we have 



In the following four figures, the theoretical cross sections from Section [2] are scaled with the 
summation factor on the right side of this equation, and plotted, along with the numerical 
data, against A (or D = S(e 2 )). Note that the summation factor has to be computed 
separately for each theoretical expression. 

Exactly the same prescription can be applied to determine the cross section for eccen- 
tricity changes. The resulting differential cross sections are given in Figs. [IT] and d2] for hard 
and soft planets, respectively. For comparison with analytic expressions, we also plot in these 
figures a number of theoretical results from Section [2T2l normalized as in equation (|33|) . 




(33) 
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Fig. 11. — Differential cross section for eccentricity changes of hard planets, as a function 
of eccentricity change (strictly, D = S(e 2 )) for the hybrid Monte Carlo model. The analytic 
results obtained by formulae from Section 2 are plotted as straight lines, while the binned 
numerical results appear as points; see detailed explanation in the main text. 
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Fig. 12. — Same as Fig. [TU but for soft planetary systems. 
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Fig. 13. — Scaled differential cross section for relative energy changes of hard planets, as 
a function of relative energy change for hybrid Monte Carlo model. The analytic results 
obtained by formulae from Section 2 are plotted as lines, while the binned numerical results 
appear as points, see detailed explanation in the main text. 
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Fig. 14. — Same as Fig. [131 but for soft planetary systems. 



4-3.3. Differential cross section for eccentricity evolution 

We first consider the case of hard planets (Fig JTTil . The two theoretical cross sections in 
the figures correspond to the limits of parabolic (dotted) and extremely hyperbolic (dashed) 
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encounters. The former correspond to closer encounters and therefore larger values of \D\. 

V 2 a\ 3/2 

The crossover between the two formulae should occur at D ~ 0.3 I — — . The trend of 



Gm J 

points with \D\ ~ 10~ 4 indeed steepens appropriately as \D\ decreases below about 0.01. 
However, since there is a considerable range in the semi-major axis, this figure merely in- 
dicates that the range over which the change of slope takes place is reasonable. Figure [JJ 
shows that there is a large overlap (in r v j a) of parabolic and hyperbolic encounters, and so 
the position of the data points between the two theoretical lines is acceptable. The notable 
exception is the strong flattening below \D\ ~ 10~ 4 . This saturation results, however, from 
the chosen value of p ma x- 

Next, we turn to the case of soft planets (Fig. [T2]) . where the interpretation, in terms 
of tidal impulsive encounters, is more straightforward. The plotted analytical result (equa- 
tion (fTTj) ) is satisfactory down to values below D ~ 10~ 4 . The saturation value in the small 
D limit is entirely consistent with that expected, using equation for encounters at a 
distance r p ~ p max . 



4-3.4- Differential cross section for fractional energy change 

Differential cross sections for the relative energy changes are displayed in Figs. [13] and 
[TH Here it is the case of hard planetary systems which is most straightforward. Almost all 
of the data agrees very well with the analytic expression for hyperbolic adiabatic encoun- 
ters (equation (1317]) ). The apparent cutoff at |A| = 1 is spurious, since the validity of the 
formula is restricted to large values of |ln|A|| (or |A| << 1. Figure [1] indicates that the 
closest encounters should be nearly parabolic, and so we have also plotted the correspond- 
ing theoretical result (dashed in Fig. [T3]) . It is not clear whether it plays a significant role 
here. We do notice, however, that softening encounters (with A < 0) have a higher cross 
section than hardening encounters, in accordance with the qualitative discussion at the end 
of Section I2XT1 

We attribute the rise in the numerical cross section at small | A| in Figure [TBI to numerical 
errors. For encounters with p ~ p ma x, we estimate from equations (12BT) and that A ~ 
10~ 44 (where even the exponent is an order-of-magnitude estimate!). This magnitude is much 
smaller than numerical errors, and even much smaller than the start-up error of initiating 
the three-body integration at a finite distance of order p ma x- Such encounters thus have 
errors which place them at the left hand side of this figure. 

The interpretation of the cross section for energy changes in soft systems (Fig. UM is 
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the most complicated of all. It is also most relevant since dissolution occurs more frequently 
in this limit. At the extreme left of this figure, the result of numerical errors is again visible. 
Over a range of small values of |A|, the numerical results lie just below the theoretical 
curve for adiabatic hyperbolic encounters (equation ([3"Uj) ). This |A| dependence intersects a 
steeper curve which represents the impulsive tidal encounters (equation ({18]) ). The numerical 
results follow the steepening trend. By comparing the two equations we expect the crossover 

2tt /Gm^ 3/2 



to occur at lAlnlAII ~ — — - . Though most of the encounters with "soft" planets 

3 \aV 2 J 

give values of V 2 a/(Gm) above 1, the median value is only of order 3 (see Fig. [2]), and so 
the crossover is expected to occur at values of |A| no less than about 0.1. It seems likely 
that the interpretation of the data is complicated by a considerable admixture of encounters 
which really should be classified as hard. 

For relatively large values of |A| ~ 0.1, we differentiate positive and negative energy 
changes (see equation (I22l) ) by multiplying the value given in equation (fl8|) with a factor 
(1 + 3 1 A|) (Fig. UM . In magnitude this corresponds quite well to the separation in the data 
points between those for A > (circles) from those for A < (squares). 

Also shown for soft planets are the two differential cross sections for impulsive non- 
tidal encounters (equations ff23l ). Though they again display the difference between positive 
and negative values of A, it is not clear what quantitative role they play in relation to 
this data. We presume that this is due to the fact that the systems are not really very 
soft. By comparing equations (123|) and ({TBI) , we estimate that the crossover should occur at 

A ~ -=— . It is only for considerably softer systems that one would expect to see this 

Vy/a 
crossover. 



5. Summary and Discussion 

5.1. Methodological Aspects 

Very few planets are found in old open and globular clusters; we have explored the 
scenario that this deficiency is associated with disruption of planetary systems in long-lived 
star clusters due to stellar encounters. Freely floating planets observed in star forming 
regions support our scenario. We have studied two prototype systems with direct A-body 
(NBODY) and hybrid Monte Carlo (HMC) simulations of star clusters which contain a large 
number of planetary systems. The direct A-body model represents a cluster similar to the 
Hyades with 20,000 stars and 1000 planets, while the HMC models represent a centrally 
concentrated globular cluster such as 47 Tuc with 300,000 stars and 30,000 planets. We 
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have recorded changes of planetary orbits due to encounters and the creation of free floaters 
(remaining in the cluster) and escapers. Our main results are summarized in Tables El [3j 
We have shown that there is a consistency between analytic and numerical models in terms 
of the cross section of encounters between planetary systems and a single star, and, to first 
order, the two models give consistent results for the liberation rates of planets (creation of free 
floaters). Encounters have been measured and recorded in large numbers and compared with 
analytic estimates for the changes of semi-major axis and eccentricity. It turns out that the 
majority of encounters in our simulations are well approximated by the analytic results, which 
corresponds to the fact that the encounters are either adiabatic (hyperbolic or parabolic) or 
impulsive and tidal in nature. In the adiabatic regimes, evolution of the planetary orbital 
elements is a diffusive process and proceeds in both directions, i.e. with similar probabilities 
in eccentricity excitation and damping as well as semi major axis increases and declines. 
However, there is also a non-negligible number of encounters which are not adiabatic and 
lead to relatively small minimum distances. These close encounters give stronger changes 
in the orbital parameters of planetary systems, and in particular a net injection of orbital 
energy and increase in the semi major axis. 

Two fundamentally different numerical methods are used in this series of investigations. 
They are the direct iV-body and Hybrid Monte Carlo algorithms, the latter being based 
on the Fokker-Planck equation for the single stars. Both schemes provide consistent results 
with each other and the analytic results (see figures and Tables). Remaining differences 
e.g. in the rates of liberation of planets per crossing time, are attributable to the more 
complex dynamical nature of the iV-body system - for many encounters with larger impact 
parameters the assumption that they are uncorrelated isolated two-body encounters is no 
longer fully justified. The examination of how well diffusion coefficients for individual two- 
body encounters are reproduced in direct iV-body systems is, howe ver, beyond th e scope of 



this paper. Interested readers are referred to discussions elsewhere ( iTheund 119961 ) . 



5.2. Application to real clusters or associations 

In order to allow the reader to extract more specific information from our general results, 
we are presenting some data regarding the free floaters in our simulations, and a scaling 
procedure to estimate the corresponding results for some real clusters. 

In the following Table H] we provide the timescale ts for a planet to become a free floater, 
by using data of Tabled averaged between the two physically equivalent runs). We obtain 
Tg by the definition rg = r cr /x p i„Ff,cr (or Tg = r r h/a;pi-fr,rh) which should give an equivalent 
result except for roundoff errors). Table H] provides exemplary values for tr/t ct and Tff/r r h, 
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Table 4: Time Scale Tg for the liberation of free floaters relative to r cr and r r h 



Model 


W r cr 


Ts/Trh 


t s Gyrs (Hya) 


r ff Gyrs (ONC) 


Tff Gyrs (47 Tuc) 


HMC soft 1 


1128.0 


1.239 


15.79 


0.367 


0.92 


HMC hard 1 


54054. 


59.52 


756.8 


17.57 


44.27 



respectively, which in fact are just the inverse averaged x-values taken from Table [2j 

To determine values of Tg in physical units we must employ the scaling of Tg with cluster 
parameters. If we assume planetary semi-major axes to stay constant in model units, as was 
done in the scaling between HMC and iV-body models here, Tg is independent of N, just as 
£pi-ff,cr- So if a typical half-mass crossing time can be deduced from observational data for 
some real clusters (as we will do below for the Orion Nebula and the Hyades clusters) we 
can just directly compute Tg from that and our measured value of Xpi_ff iC r- If, however, a 
value of the half-mass relaxation time is more well known from observations (as in the case 
of the globular cluster 47 Tuc), we have to be aware that x p i_ff ir h scales with the relaxation 
time (see explanations after Table [2]). 

The data to obtain lifetimes of planetary systems in physical units are obtained as 
follows: The globular cluster 47 Tuc is a typical old stellar cluster with a very high central 
density (of order 10 6 M Q /pc 3 ). We expect this cluster to be the most hostile environment for 



planetary systems. From the Harris catalog (IHarria Il996l ) we get the following data for 47 



Tuc: Th = 2.79', with a distance of R = 4.5 kpc, which implies that =3.6 pc; the half-mass 
relaxation time is given as logt r h = 9.48, i.e. = 3.02 Gyrs. In order to apply our model, 
we need to determine the particle number, or total mass, which are poorly constrained in 
the literature. Therefore we use the following procedure, which is sufficient for an order-of- 
magnitude approach. We adopt the relaxation time scaling of t r h oc Nr^ / log^N). From 
our HMC model, with N = 300.000, r h = 0.77pc, and t lh = 7.28 • 10 7 yrs, we then deduce 
a particle number of N = 1.4 • 10 6 stars for 47 Tuc. The scaling factor for the relaxation 
time between our HMC model cluster and 47 Tuc is 4.06. Hence we predict a time scale Tg 
to be 0.305 (14.66) relaxation times for a planet to become a free floater, i.e. 0.92 (44.27) 
Gyrs. These results have been obtained from our finding that Tg / t r h scales inversely with the 
relaxation time, because in systems with larger N more planets are liberated per relaxation 
time). The values outside and inside brackets are for the soft and hard planetary systems 
respectively. In other words, we predict for 47 Tuc, the dissolution rate is all (7 %) of soft 
(hard) planetary systems per relaxation time. We have said that the hard time scale is 14.66 
relaxation times (44.27 Gyrs, based on Table Hj); here "soft" refers to planets with semi-major 
axes between 3 and 50 AU, while " hard" refers to 0.03 to 5 AU, both in logarithmically equal 



distribution. iFregeau et al.l (120061 ) quote a lifetime of 10 Gyrs for hard planetary systems 
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(0.05 A U) in their statisti c al an alysis, which is also in qualitative agreement with earlier 



work of ISmith &: Bonnell I (120011 ). From these previous papers, one would have expected 
that most present-day planetary systems in globular clusters have been dissolved. While the 
given t ime scales are o nly a rough estimate, detailed models as presented in this paper and 
also in iFregeau et al.l (120061 ) suggest that a significant fraction of planets wi th semi-major 
axes b etween 0.03 and 5 AU (hard in our notation, intermediate in that of IFregeau et al. 



( 120061 ) ) could have survived stellar encounters in this cluster (see further discussions below). 



The Orion nebula cluster (ONC) is an example of a young star cluster with intermediate 
central density, of order 2 ■ lO 4 M /pc 3 . For this system, we use a half-mass radius of 0.8 pc 
and a velocity dispersion of 2.34 km/s, to obtain the half-mass crossing time as 3.25 ■ 10 5 
yrs. For ONC also the h alf-mass relaxation time (t r h = 6.5 • 10 6 yrs) and particle number 
N = 2800 are provided ([Hillenbrand fc Hartmannl 119981 ) . Planetary system lifetimes are 



then obtained by using our measured lifetimes in units of crossing times (see Table Hj). 

Note that a scaling using the measured half-mass relaxation time of the ONC would 
provide approximately the same results, if the proper Coulomb scaling factor is applied to 
our measured time scales (in units of relaxation time, see Tables [2J [3]) . 

For the Hyades cluster, which is a relatively low density open cluster we use the half- 
mass radius = 4.5 pc and central velo city dispersion of 0.3 km /s to determine a half-mass 
crossing time of approx. 1.4 ■ 10 7 yrs (jPerryman et al.lll998l ). Again we derive physical 
timescales in Table SJ Note that these estimates are subject to many uncertainties, such as 
the use in our models of stars of equal mass. 



5.3. Total Cross Sections and Timescales to get Free Floaters 



Fregeau et al.l (120061 ) provide in their Monte Carlo simulations cross sections for orbital 
changes of planetary systems obtained from a large set of simulations of three-body en- 
counters. We want to compare their results with ours for the case of ionization of planetary 
systems due to relatively close encounters. For that we use the usual ansatz for the planetary 
system lifetime as 

r ff = -^-tt (34) 
na ion V 

where n = p/m is the average particle density, a lon is the total cross section for dissolution 
of planetary systems, and V is the typical encounter relative velocity. The relevant total 
cross section for ionization in our case is taken from equation ( |23~1) for non-tidal impulsive 
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encounters as 

87iGm r „,. N . N1 
o"ion = y 2 a [F(A max ) - F(A min )J 

f < A > = K^" 1 ) (35) 

where we have obtained the total cross section for "ionization" by integrating equation (1231) 
(A < 0) from some suitable A min to A max . With A max — ► — oo and the assumption that 
A m i n = — 1 (the minimum required relative energy change for ionization is unity) we get 
F(A max ) = 11/3, thus 

40vr Gm 

o"ion = -g — y^a (36) 

and 

Tff = 40^ (37) 



This result is in agreement with iFregeau et al.l (120061 ) in the case of equal stellar masses 



although we get for the case of 47 Tuc slightly larger values of the planetary lifetimes. 
Possible reasons are that the distribution of encounter velocities in a real cluster differs from 
those in three-body experiments, and that planets, which just escape from a three-body 
encounter with a very small velocity at infinity could easily be recaptured by their own host 
star or another star due to potential fluctuations. 

Note that the lifetime of planetary systems as given in equation ( 1371) does not vary with 
the particle number N, provided the planetary system's semi-major axes remain constant 
relative to the iV-body scale (system scaling radius), consistent with our numerical results. 
If one is interested in ionisation time scales relative to the relaxation time we have to look 
at the quantity 

1 T ff GmlogA 

= jj- (38) 

£ P i-ff,rh r vh V z a 

where log A is the Coulomb logarithm contained in the relaxati on time. If we take th e 
standard Coulomb logarithm logA = log(7iV), with 7 = 0.11 (cf. Giersz fc Heggie ( 19941 ) ) 
we get for the scaling 

_7rh V 2 a N 

Xpl ~ S ' rh ~ rff ~ Gmlog( 7 iV) ~ log( 7 A0 1 J 

This is indeed the empirically measured scaling behaviour, as we have commented in the 
caption to Table [2j Note that in these considerations, also in the choice of initial conditions 
for our runs with different particle number (HMC and iV-body) we have neglected the possible 
influence of a separate variation of the planetary orbit scale a with respect to the iV-body 
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system's scaling radius R. Since rg ~ R/a the effect would be easy to estimate to first order 
(if only the size of planetary orbits relative to the cluster size is varied). If however, the scale 
changes are due to different cluster models (e.g. King models of different concentration) 
further numerical studies would be required. 

5.4. Discussion 

We have provided a further step towards a self-consistent modelling of the origin of 
differences in the frequency of planetary systems between star clusters and the galactic 
field due to encounters of planetary systems with passing stars. For two representative 
star clusters with particle numbers resembling open and globular clusters, respectively, we 
have tested analytical cross sections against two different series of self-consistent numerical 
simulations under realistic conditions of a star cluster (not merely a series of three-body 
scattering experiments). For the illustration of our results, we have provided liberation 
time scales for planets (to become a free floater) in their respective cluster environment. In 
contrast to previous work we have extended our study to planetary systems ranging from 0.05 
AU to 50 AU in semi-major axis and all eccentricities. To ease discussion and understanding 
we divide our planetary systems in two classes (hard and soft), and determine the average 
time scale for their survival, see Table HI We find that the liberation rate of planets per 
crossing time is constant, and per relaxation time scales with the scaling factor iV/logA, 
where log A is a Coulomb logarithm (provided the planetary system's size remains constant 
relative to the star cluster radial scale). 

The results have been applied to three typical cases of clusters, the high-density old 
globular cluster 47 Tuc, the intermediate density young star forming cluster in the Orion 
nebula (ONC), and the low-density open cluster of the Hyades. We find that even in the 
dense environment of 47 Tuc about one quarter of the short-period planets survive even 
after 10 Gyrs. Our results are interesting for a number of reasons. First of all, the diffusive 
nature of encounters between planetary systems and stars has been confirmed. With the 
exception of the few close encounters most intermediate and distant encounters have no 
preferred sign of change for the orbital parameters of the planet. Therefore planets can be 
scattered onto a more eccentric orbit, even if they are initially very close to the host star. 
Subsequent evolution of such planets by tidal effects may lead to their inflation and tidal 
disruption. This process is expected to be particularly effective in dense globular clusters. 
Therefore it is not guaranteed that the fraction of surviving systems we have found for the 
case of 47 Tuc would be observable - it may just be destroyed by secondary processes such as 
tidal interactions. Here we have just rescaled from our exemplary model simulations to the 
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particle number and mass of 47 Tuc; a detailed case study of 47 Tuc would require a model 
using the half-mass radius and concentration parameter properly tailored so as to reproduce 
the presently observed values. Also, a more detailed case study of 47 Tuc would require 
as well to include a stellar mass spectrum, stellar evolution, some binary fraction. This is 
subject of future work and beyond the scope of this first exploratory study. 

In contrast to the less-rich Taurus star forming region, the ONC contains massive stars 
similar to the hypothetical progenitors which generated the radioactive isotopes prior to the 
formation of the solar system. This and other, slightly less rich associations are also likely 
to disperse in ~ 100 Myr rather than evolve into an open cluster. The main issue here is 
whether the planet-bearing stars can preserve their companions before joining the field star 
population. Our numerical results indicate that in the ONC a few times 10 8 years is a critical 
threshold for the survival of wider planets and planetary systems such as the solar system. 
Thus, the disruption of some planetary systems is expected, especially near the region where 
dense protostellar cores are concentrated. These are also the regions where the existence of 
free floaters have been reported. Our models suggest that a significant fraction of the free- 
floater population may represent the relics of disrupted planetary systems. Our results also 
indicate that the velocity dispersion of the freely floating planets must be generally small, 
as most are retained by the cluster potential, at least for the case of soft planets (Tables [2] 
and ED. 

Finally, for open clusters such as the Hyades, there is no problem in preventing all the 
planetary systems (wide or compact) from destruction before the cluster dissolves on a Gyr 
time scale. The observed dichotomy in the fraction of stars with planets between the field 
stars and the Hyades requires additional explanation. The present-day age of the Hyades is 
THya — 0.8 Gyr. It is likely that the cluster was somewhat denser in the past which would 
shorten the planetary-system disruption time scale below that listed in Table HI A more com- 
pact young Hyades would also reduce its r ff below its lifespan. Both eccentricity excitation 
and modest migration would shake up otherwise stable planetary systems. If systems are 
driven into nearly resonant configurations, modest eccentricity excitation can render them 
dynamically unstable. We have not yet taken into account the dynamics of multi-planetary 
systems in our model, but this subject will be addressed in our future investigations. 

Last, but not least, stellar encounters with planetary system can alter the eccentricity 
of orbits with relatively large semi-major axes efficiently, and also they are the only process 
which can create significant amounts of inclinations out of the original orbital plane of the 
planetary system (though we have not considered that issue in this paper). These properties 
are potentially important for understanding the formation of Kuiper-Belt objects. 

We have neglected for the present paper the effect of stellar binaries and multiples, as 
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well as that of a stellar mass spectrum. In particular the presence of a large number of binaries 
would be expected to have a profound effect through the action of binary system - planetary 
system encounters. The inclusion of binary stars would certainly reduce the magnitude of tq 
and the key issue is by what extent. Also a further study of the free floaters in globular and 
open clusters is interesting. Due to their very small mass they cannot receive large amounts 
of energy in the mass segregation - equipartition processes structuring the cluster during 
its dynamical evolution. Therefore large numbers of free floating planets should be retained 
in the cluster, a conclusion in conflict at least with present observations. From the point 
of view of escape, free-floating planets would behave very much like low-mass stars. There 
is much observational and theoretical evidence that low-mass objects escape preferentially 
from star clusters, especially in the presence of steady and time-dependent tidal fields. Mass 
segregation tends to drive such objects to large radii, where they may be efficiently removed 
by tidal effects. To study the fate of a component of very low-mass planetary objects during 
the standard evolution of star cluster is another interesting future project. 

We conclude that, for an understanding of the diversity of planetary systems, the fact 
that they originate in a dense star cluster cannot be neglected. Our solar system could have 
formed in an open cluster or in the outskirts of an object which evolved from an ONC like 
star cluster. Other processes shaping planetary systems, such as resonances, and interaction 
with gas, need to be distinguished from the effect of diffusive encounters. 
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A. Formulae of three-body scattering for planetary problems 

We present formulae for the change in binding energy and eccentricity of a binary as a 
result of scattering by a distant third body. This is a comprehensive collection of main results 
from earlier papers (Heggie 1975, Heggie & Rasio 1996, Roy & Haddow 2003, Heggie 2005), 
complemented by a few new results, in order to match the requirements for comparison with 
our numerical simulations. 

These formulae assume that the encounter is adiabatic, but take account of the hyper- 
bolic geometry of the relative orbit. In adiabatic encounters we assume that the encounter 
between the host and field stars occurs on a much longer time scale than the planet's or- 
bital period. Modification of the perturbing potential due to the passage of the field star is 
essentially adiabatic. Consequently, changes of the eccentricity 5e can be computed with an 
orbit-averaging method, in which incremental change of 5e per orbit can be evaluated for 
an instantaneous field star position. For A, however, the analytical technique is a bit more 
subtle (see the original papers). 



A.l. Change of eccentricity 

A. 1.1. Non- circular binaries 

Long ago Heggie (1975), equation (5.66), derived a formula for the change in eccentricity 
of a binary subject to a parabolic flyby of a third body. Corrected for an overall sign error, 
it is: ^ 2 

5e = - 15l L - ^ eVT^(-) (a-Ab-A + a-Bb-B), (Al) 



4^2 VM 12 M 123 \r P J 



This preprint was prepared with the AAS IATeX macros v5.0. 
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where e is the eccentricity, is the mass of the perturber, Mi 2 = vri\ + m 2 is the total 
mass of the binary, M 123 = vn,\ + m 2 + is the total mass of the system, a is the (initial) 
semi-major axis of the binary, r p is the distance of closest approach between the perturber 
and the centre of mass of the binary (on a Keplerian approximation), a is a unit vector along 
the pericentre of the binary, b is an orthogonal unit vector in the plane of motion of the 
binary (so that a A b i is directed along the angular momentum vector of the binary), A is a 
unit vector along the pericentre of the third body, and B is an orthogonal unit vector in the 
plane of motion of the third body (so that A A B is directed along the angular momentum 
vector of the relative motion of the third body and the binary). Heggie & Rasio (1996) gave 
the corresponding formula for a hyperbolic flyby: 



5e 



5e 



m|a 3 (l 



x 



2r " \/ M 12 M 123 r3(l + e 
a Ab • A 



^3 



3e' arccos ( 

e 



+a - Bb B 



■j i 2 I 1 

3e arccos 



+ (4e /2 - l)Ve' 2 - 1 
(2e /2 + l)v / e' 2 - 



(A2) 



where e' is the eccentricity of the third body. Expressing the unit vectors in the equation 
above using orbital elements as in Roy & Haddow (2003), equation (18) (angles u, Q, i) we 
can write the result as follows: 



8e 



15, 



rrin 



i/2/a\3/2 ey /i 
r 



(1 + e'fl 2 



x 



4 VM 12 M 123 
x | sin 2 i sm(2Q) arccos(— 1/e') + Ve' 2 — 1 

+ - [(1 + cos 2 i) cos(2u;) sin(2fi) + 2 cosi sin(2u;) cos(2Q)] 
3 



.-/2 



1) 



3/2 



(A3) 



where Q is the longitude of the ascending node of the orbit of the third body, measured in 
the plane of motion of the binary from a particular origin; i is the inclination of the two 
orbits; and u is the longitude of pericentre of the third body, measured from the ascending 
node. Thus we have arrived at equation (7) of Heggie & Rasio (1996). In Section 2.2 and 
Appendix [B] we use this expression to estimate cross sections. 



A. 1.2. Circular binaries 



The above formulae are all that is needed for most purposes, but for the sake of expo- 
sition we also give here two formulae, also from Heggie & Rasio (1996), for the case when 
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e = initially. For hyperbolic encounters the result is: 



Se = 3v^^^ 3/4 < e ' + .|» 3/4 x 



M 5/4 

JW 123 



x exp 



x cos 



M 



12 



1/2 



M 123 
, i 4 



r p \3/2 



arccos(l/e') 
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(e' - l) 3 / 2 
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% 4 2 2 
cos" — I — sin 4 — I — cos 2 - sin 2 - cos(4u> + 2Q) 
2 9 2 3 2 2 v ' 



1/2 



This result simplifies a little for parabolic motion of the third body, to 

l / 2 ,^ n 3/2 



<5e 



3V2tt 



M 



5/4 
123 



exp 



2 / 2M 12 

3 V M 123 



X 



X COS 



4-4*4 2 1 ■ 2 1 



(•( >s" 77 + 77 sin" 1 75 + 77 cos - sin - cos(4w + 2J7) 



1/2 



(A4) 



(A5) 



A. 2. Change of binding energy 

A. 2.1. Non- circular binaries 

Roy & Haddow (2003) give the following expression for the change in energy of the 
binary in the case of a parabolic flyby: 



Se = - Gmim2m3 ^^ 5/2 exp(-2K/3)x 



M 12 r3 120" 



+ 



x |60a 2 ei (sin(nt )(a.B) 2 + 2cos(nt )a A a B - sin(nt )(a • A) 2 ^ 

+ 120a6e 4 (— sin(nt )b • B a • A — sin(nt )a • B b A+ 
+ cos(nt )a • B b • B — cos(nt )a ■ A b ■ A) + 

+606 2 e 2 (sin(nio) (b • A) 2 - 2 cos(nt )b B b A - sm(nt ) (b • B) 2 ) | , (A6) 



where K 



'2M 12 r 3 





M 123 a 3 



a\/l — e 2 is the semi-minor axis of the binary, n 



GM l2 



is 



the mean motion of the binary, to is the time of pericentric passage of the binary (on a time 
scale in which closest approach of the third body occurs at t — 0), and the coefficients are 
given by: 



ei = J-i(e) - 2eJ (e) + 2eJ 2 (e) - J 3 (e) 
e 2 = J-i(e) - J 3 (e) 
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e 3 = eJ_i(e) - 2 J (e) + 2J 2 (e) - eJ 3 (e) 

e 4 = j_ 1 ( e )_ e j ( e )_ e j 2 (e) + J 3 (e); (A7) 
here J n is the Bessel function of the first kind of order n. 

By combining the mathematical procedures in the two papers Heggie & Rasio (1996) and 
Roy & Haddow (2003) it is not hard to show that the corresponding result for a hyperbolic 
encounter of eccentricity e' gives the following expression: 



5e 



Gmim 2 m 3 ^phx ( M 12 \ 5/4 ^r p \ 3/4 ( e ' + l) 3 / 4 



Mi 2 a 3 120 V. M i23/ e' 2 
exp ^ e' 2 — 1 — arccos — J ^ |60a 2 ex ^sin(72t )(a ■ B) 2 + 

+2 cos(nt )a • A a • B — sin(nt ) (a • A) 2 J + 

+120a6e 4 (— sin(nt )b • B a • A — sin(nt )a • B b • A+ 
+ cos(nto)a • B b • B— 

- cos(nt )a • A b • A) + 60b 2 e 2 (sin(nt )(b ■ A) 2 - 

-2 cos(nt )b B b A - sin(nt ) (b • B) 2 j | ; (At 



here n' = \l — — where a' is the semi-major axis of the hyperbolic motion, and so 



a' 3 



r p = a!(e' — 1). Now expressing again the result by using orbital elements we get: 

/ e' 2 — 1 — arccos(l / e') 



G mi m 2 m 3 ^ (e' + l) 3 / 4 5/2 / K ( V~e 



!)3/2 



X 



p 

x |eia 2 [sin(2u; + nto)(cos(2z) — 1) — sin(2co> + nto) cos(2f2) cos(2i) — 

—3 sin(2u; + nto) cos(2fi) — 4 cos(2a; + nto) sin(2f2) cos i] + 

+e 2 6 2 [sin(2a; + nt )(l — cos(2i)) — sin(2a; + nto) cos(2f2) cos(2i) — 

—3 sin(2co> + nto) cos(2fi) — 4 cos(2u; + nto) sin(2f2) cos i] + 

-\-e±ab [—2 cos(2c<j + nto) sin(2fi) cos(2i) — 6 cos(2a; + nto) sin(2fi) — 

-8sin(2u; + rato)cos(2n)cos2]} (A9) 

where Q is measured from a. From here the interested reader should read off the definition 
of the functions F, fi, and f 2 used in equations (T2"4"l) and ([25]) . 



A. 2.2. Circular binaries, parabolic case 



Roy and Haddow (2003) also provided an expression for the energy change of a circular 
binary in the case of parabolic flyby This formula, expressed in terms of the orbital elements, 
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is as follows: 

Se = g r4 ^ 7/2 exp ^--ifj (/i 2 - ^)(1 + cos?) sin 2 z x 

x [(cos 3 uj — 3 sin 2 u; cos cu) sin(nt ) + (3 cos 2 u sin u — sin 3 a;) cos(nt )] (A10) 

where \ii = rrii/M^ {i = 1, 2) are the reduced masses. 



B. Cross sections for the change of eccentricity 
B.l. Impulsive tidal encounters 

Here our starting point is an expression for the total change, due to the passing star, i n 



the specific angular momentum J of the planet relative to its sun. Following iHeggid (119751 ). 



where (in Section 4.2) an analogous calculation of the energy change is given, we obtain 

Here p is a vector from the sun to the closest point of the (rectilinear) path of the passing 
star, and r^ is the component of the position vector of the planet, at this time, which is 
orthogonal to the path of the passing star. As usual we have specialised to the case of equal 
stellar mass and negligible planetary mass. 



In the tidal regime we expand in powers of rj_, obtaining to lowest order 

2Gm 

where p is the unit vector parallel to p. Hence 



5J = —v x (-r ± + 2(r_Lp)p) , (B2) 



2Gm 

J-5J = ~Yf v - r x (r x (_r± + ( r± -P)p))- (B3) 

The corresponding result for the change in energy, 5e, is similar, except that the factor 
rx in equation (1B2|) is replaced by v ■ , where v is the velocity of the planet (relative to its 
sun) at the time of closest approach of the passing star. By equation (jSJ) the change in e 2 
can now be written as 

D = s^ A -(- rl + 2r '^ < B4 » 

where A = J 2 v + 2eJ x r. Now Aj_, r_|_ and p are vectors in the same plane, and so D = 
4r_i_v4j_ 

— tt~o cos(ip + 29), where 9, tb are the angles in this plane from rj_ to p, Aj_, respectively. 
GmVp 2 
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Thus the differential cross section, before averaging all the binary parameters, is 



%=.h{ D+ ^-^ +2 °^°- < b5 > 



Hence 



da f pdp 



\GmVp 



2r ± A ± 
GmVD 2 ' 



(B6) 



(B7) 



Next we note that r± scales as a, the semi-major axis, and A scales with (Gm) 3/,2 a 1//2 . 
Hence we can write 

da 2(Gm) 1 /2 a 3/2 



dD VD 2 
the extra subscripts denoting these normalisations. 



■r±iA ± i, (B8) 



This is as far as we have been able to take the analytic evaluation of this cross section. 
Resorting to Monte Carlo integration to average over the normalised binary parameters, we 
found that 

da 2{Gm) 1 / 2 a^ 2 C l 



dD VD 2 
where C\ = 0.883 to three significant figures. 



(B9) 



B.2. Extremely hyperbolic adiabatic encounters 

We begin with equation (1A2I) . let mi = = m and rri2 = 0, and take the asymptotic 
form for e' — > oo. In this limit r p ~ p, the impact parameter, and we evaluate the limiting 
expression for e' using equation ([2]). The result is 

heJGmil - e 2 )a 3 / 2 , A j ^ j A\ 
5e = ^ ^— 2a.Ab.A + aBb.B , BIO 

Vp 2 

where a, b are certain orthogonal unit vectors in the plane of motion of the planet, and A, B 
are similarly defined in the plane of motion of the passing star (Appendix I A. 1 . 1 1 ) . Hence the 
total cross section for D is, for D > 0, 



9 l0ire 2 JGm(l-e 2 )a 3 / 2 . A 4 - s a . 
a{D > D ) = np 2 = v '- (2a.Ab.A + aBb.B), (Bll) 
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provided that the last factor is negative. 



Now we average over the orientation of the planetary orbit. To do this we think of A and 
B as providing the first two vectors of a basis, and the expression to be averaged is therefore 
X = 2a\b\ + d2&2, where this is negative. With the usual spherical polar coordinates we take 



(sin 9 cos 6, sin 9 sin d>, cos ( 



(B12) 



and then b, which is an arbitrary unit vector orthogonal to a, may be written as 

b = (cos 9 cos cos if) — sin sin if), cos 9 sin cos if) + cos sin if), — sin# cosip). (B13) 
Therefore X = a cos if) + (3 simp, where 



a = sin 9 cos 9(2 cos 2 + sin 2 - 
f3 = — cos sin sin ^. 



Averaging with respect to if) for X < gives 



(X) = — V« 2 + /3 2 - 



7T 



(B14) 
(B15) 



(B16) 



It appears that averaging over 9 and has to be done by numerical quadrature, and yields 

2C 2 



(X) 



7T 



(B17) 



where C2 = 0.5932 approximately. 



Finally we average over a thermal distribution of e, and differentiate with respect to D 
to obtain the differential cross section 



da I6C2 V Gma 3 
dD = ~3 VD 2 
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B.3. Near-parabolic encounters 



We begin with equation (19) in lHeggie &: Rasiol (1l996l ). Specialising to the case of two 
equal masses m and one vanishing mass, this becomes 



. . 9v/3 /15tt\ 2/3 
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(B19) 
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Eliminating 5e in favour of D = S(e 2 ) ~ 2e5e, and then differentiating with respect to D, 
we arrive at the differential cross section 

1 2 



da 3\^3 



dD 



7tt 



(15tt) 



2/3 



Gma 



,4/3 



(l-e 2 ) 1 / 3 ^) 



-5/3 



Finally, averaging over a thermal distribution of eccentricities, i.e. /(e) = 2e, gives 

2 



dD 21 v ; 



Gma 



-5/3 



(B20) 
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C. Cross Sections for the Change in Energy 

Here we consider only adiabatic encounters, and turn first to the parabolic case. Equa- 
tion (1A6I) shows that 

2iT 



A = gexp (—J, (CI) 

where = (r p /a) 3//2 for the present assumptions about the masses, and g is a dimensionless 
coefficient which depends on K and on the geometry of the encounter. Remarkably, we do not 



need to evaluate it in order to obtain the cross section to lowest order. Using equation (IA6I) 
in the parabolic limit, we deduce that 



V V 2 gj 

where the factor of 1/2 in the second expression is needed because the geometric factor g 
has the correct sign in only half of all encounters. Differentiating with respect to A, we can 
ignore the derivative of g; g depends on powers of r p or p and hence on powers of In A, as 
we can see from equation (IC2I) . and so it varies much more slowly than A. Similarly we can 
neglect ln|g| by comparison with |A| for sufficiently weak encounters. Therefore 

da 2irGma 1 



dA " V 2 llnlAllVsiAl 



(C3) 



The extremely hyperbolic case follows a similar line of argument. As stated in equa- 
tion ( 1271) we have 

A = gexp(-^Y (C4) 



3V^ 

Now we use equations (T5]) and (jl]) to express everything in terms of p, and find that 
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Differentiation with respect to A with the same approximations leads to the result 
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